Input

Analysis for FDCR manuscript describing the temporal dynamics of drug cue reactivity with replication. ## Read Data

#function to be able to read in an .RData file and assign it to a variable
loadRData <- function(fileName){
  #loads an RData file, and returns it
  load(fileName)
  get(ls()[ls() != "fileName"])
}
idps_neurocaps <- loadRData('../paper-dcr-temporaldynamics/data/neurocaps/idps-neurocaps-2-13-2020.RData')
idps_neurocaps$motion <- idps_neurocaps$dcr_motaveraw_r1
#remove subjects due to excessive motion
idps_neurocaps <- idps_neurocaps[idps_neurocaps$motion < 0.3 & (!idps_neurocaps$neurocaps_exclude),]
idps_neurocaps_ocr <- idps_neurocaps[idps_neurocaps$dcr_version == 'ocr',]
idps_neurocaps_mcr <- idps_neurocaps[idps_neurocaps$dcr_version == 'mcr',]
idps_tacs <- loadRData('../paper-dcr-temporaldynamics/data/tacs/idps-tacs-2-13-2020.RData')
idps_tacs$motion <- idps_tacs$dcr_tpre_motaveraw_r1
idps_tacs <- idps_tacs[idps_tacs$motion < 0.3 & (!idps_tacs$tacs_exclude),]
idps_tdcs <- loadRData('../paper-dcr-temporaldynamics/data/tdcs/idps-tdcs-2-13-2020.RData')
idps_tdcs$motion <- idps_tdcs$dcr_tpre_motaveraw_r1
idps_tdcs <- idps_tdcs[idps_tdcs$motion < 0.3 & (!idps_tdcs$tdcs_exclude),]

#only use neurocaps OCR subjects who also did tACS
idps_neurocaps_ocr <- idps_neurocaps_ocr[idps_neurocaps_ocr$id %in% idps_tacs$id,]

#only use tacs subjects with good neurocaps data
idps_tacs <- idps_tacs[idps_tacs$id %in% idps_neurocaps_ocr$id,]

#idps_tacs$dcr_tpre_response_craving_[0-7]
#idp_tdcs
#dcr_tpre_response_box_rt_0
#dcr_tpre_response_craving_rt_0
#idps_neurocaps_ocr$dcr_response_box_craving[0-7]


n_tdcs <- length(unique(idps_tdcs$id))
n_neurocaps_mcr <- length(unique(idps_neurocaps_mcr$id))
n_neurocaps_ocr <- length(unique(idps_neurocaps_ocr$id))
n_tacs <- length(unique(idps_tacs$id))


print('Subjects In Each Dataset')
## [1] "Subjects In Each Dataset"
print('tDCS MCR')
## [1] "tDCS MCR"
print(n_tdcs)
## [1] 65
print('Neurocaps MCR')
## [1] "Neurocaps MCR"
print(n_neurocaps_mcr)
## [1] 29
print('Neurocaps OCR')
## [1] "Neurocaps OCR"
print(n_neurocaps_ocr)
## [1] 22
print('tACS OCR')
## [1] "tACS OCR"
print(n_tacs)
## [1] 22
#write.csv(idps_tdcs$id, 'subjects-tdcs.csv', row.names = FALSE)
#write.csv(idps_neurocaps_mcr$id, 'subjects-neurocapsMCR.csv', row.names = FALSE)
#write.csv(idps_neurocaps_ocr$id, 'subjects-neurocapsOCR.csv', row.names = FALSE)
#write.csv(idps_tacs$id, 'subjects-tacs.csv', row.names = FALSE)


#write out behavioral data to share with Hamed/Henry

behavioral_cols <- grepl('response', names(idps_tdcs))
behavioral_cols[names(idps_tdcs) == 'id'] <- TRUE
#write.csv(idps_tdcs[, behavioral_cols], 'behavior-tdcs.csv', row.names = FALSE)

behavioral_cols <- grepl('tpre_response', names(idps_tacs))
behavioral_cols[names(idps_tacs) == 'id'] <- TRUE
#write.csv(idps_tacs[, behavioral_cols], 'behavior-tacs.csv', row.names = FALSE)

behavioral_cols <- grepl('response', names(idps_neurocaps_mcr))
behavioral_cols[names(idps_neurocaps_mcr) == 'id'] <- TRUE
#write.csv(idps_neurocaps_mcr[, behavioral_cols], 'behavior-neurocapsMCR.csv', row.names = FALSE)

behavioral_cols <- grepl('response', names(idps_neurocaps_ocr))
behavioral_cols[names(idps_neurocaps_ocr) == 'id'] <- TRUE
#write.csv(idps_neurocaps_ocr[, behavioral_cols], 'behavior-neurocapsOCR.csv', row.names = FALSE)

#idps_tacs <- read.csv('idps-tacs-abstract.csv')
#variable names are:
#for tACS
#dcr_tpre_stats_tdcsprelim_[condition].[run].0.coef_mean_[roi]
#for neurocaps
#dcr_stats_tdcsprelim_[condition].[run].0.coef_mean_[roi]

#conditions are: drug or neutral
#runs are (e.g. run 1 block 1):
#r11, r12, r13, r14

#ROIs are:
#1: VMPFC extracted at p < 0.001
#3: LSTG extracted at p < 0.001
#4: RSTG extracted at p < 0.001
#5: L Ventral Striatum extracted at p < 0.05 intersected with brainnetome ROIs 219 and 223
#6: R Ventral Striatum extracted at p < 0.05 intersected with brainnetome ROIs 220 and 224
#8: R Amygdala extracted at p < 0.05 intersected with brainnetome ROIs 212 and 214
library(reshape2)
library(lme4)
## Loading required package: Matrix
library(lmerTest)
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
library(ggplot2)
library(sjstats) #for icc of mixed effects models
plot_one <- function(roi, this_label, idps, limits_spag, limits_y = c(-1, 1), prefix = 'dcr_', task = ''){
  #roi <- '1'
  #this_label <- 'VMPFC'
  
  
  #roi <- '1'
  #this_label <- 'VMPFC'
  #prefix <- 'dcr_tpre_'
  #idps <- idps_tdcs
  
  column_prefixes <- c('stats_tdcsprelim_drug.r11.0.coef_mean_',
             'stats_tdcsprelim_drug.r12.0.coef_mean_',
             'stats_tdcsprelim_drug.r13.0.coef_mean_',
             'stats_tdcsprelim_drug.r14.0.coef_mean_',
             'stats_tdcsprelim_neutral.r11.0.coef_mean_',
             'stats_tdcsprelim_neutral.r12.0.coef_mean_',
             'stats_tdcsprelim_neutral.r13.0.coef_mean_',
             'stats_tdcsprelim_neutral.r14.0.coef_mean_')
  this_roi <- c('id', 'motion', paste0(prefix, column_prefixes, roi))
  one_dataset <- idps[, this_roi]
  long_data <- melt(one_dataset, id.vars = c('id', 'motion'))
  long_data$condition <- NA
  long_data$condition[grepl('neutral', long_data$variable)] <- 'neutral'
  long_data$condition[grepl('drug', long_data$variable)] <- 'drug'
  #put neutral in the intercept
  long_data$condition <- factor(long_data$condition, levels = c('neutral', 'drug'))
  long_data$time <- NA
  long_data$time[grepl('r11', long_data$variable)] <- 1
  long_data$time[grepl('r12', long_data$variable)] <- 2
  long_data$time[grepl('r13', long_data$variable)] <- 3
  long_data$time[grepl('r14', long_data$variable)] <- 4

  #mean center on time
  long_data$time <- long_data$time - mean(long_data$time)


  #this_lme <- lmer(paste('value ~ condition * time + (1|id)'), data = long_data)
  #for checking for NA's, looks like it's all good now
  #print(long_data)
  this_lme <- lmer(paste('value ~ condition * time + motion + (1|id/condition)'), data = long_data)
  #print(summary(this_lme))
  within_visit_iccs <- sjstats::icc(this_lme)

  n_subjects <- length(unique(long_data$id))
  means <- aggregate(long_data$value, by = list(long_data$condition, long_data$time), FUN = mean, na.rm = TRUE)
  names(means)[1:2] <- c('condition', 'time')
  sds <- aggregate(long_data$value, by = list(long_data$condition, long_data$time), FUN = sd, na.rm = TRUE)
  names(sds)[1:2] <- c('condition', 'time')


  #print(this_label)
  mean_str <- paste(this_label, '_mean', sep = '')
  names(means)[3] <- mean_str 
  

  sd_str <- paste(this_label, '_sd', sep = '')
  se_str <- paste(this_label, '_se', sep = '')
  names(sds)[3] <- sd_str
  
  plot_frame <- merge(means, sds)
  plot_frame[, se_str] <- plot_frame[, sd_str] / sqrt(n_subjects)
  plot_frame$bar_top <- plot_frame[, mean_str] + plot_frame[, se_str]
  plot_frame$bar_bottom <- plot_frame[, mean_str] - plot_frame[, se_str]
  
  plot_frame$label <- this_label
  
  ebsize <- 1.5
  linesize <- 1
  #swap condition levels back for this plot, so that drug is red and neutral is blue
  plot_frame$condition <- factor(plot_frame$condition, levels = c('drug', 'neutral'))
  p <- ggplot(plot_frame) + geom_line(aes_string(x = 'time', y = mean_str, color = 'condition'), size = linesize) + geom_errorbar(aes_string(x = 'time', ymax = 'bar_top',                                                       ymin = 'bar_bottom', color = 'condition', width = 0.1), size = ebsize) +
                      theme(text = element_text(size=10)) + ggtitle(task) + ylab('') + xlab('') + ylim(limits_y)
  
  #print(p)
  
  
  p2 <- ggplot(data = long_data[long_data$condition == 'drug',], aes_string(x = 'time', y = 'value', group = 'id')) + 
    geom_line() + stat_summary(aes(group = 1), geom = "point", fun.y = mean, shape = 17, size = 3) + 
    stat_smooth(aes(group = 1)) + labs(
    x = "Time",
    y = "% Signal Change",
    title = "drug") 
  #print(p2)
  p3 <- ggplot(data = long_data[long_data$condition == 'neutral',], aes_string(x = 'time', y = 'value', group = 'id')) + 
    geom_line() + stat_summary(aes(group = 1), geom = "point", fun.y = mean, shape = 17, size = 3) + 
    stat_smooth(aes(group = 1)) + labs(
    x = "Time",
    y = "% Signal Change",
    title = "neutral") 
  #print(p3) 
  
  wide_icc_data <- dcast(long_data, id ~ variable)
  within_iccs_simple <- ICC(wide_icc_data[, 2:ncol(wide_icc_data)])
  within_iccs_simple_intervals <- extract_icc_intervals(within_iccs_simple)
  
  
  wide_icc_data_contrasts <- dcast(long_data, id + time ~ condition)
  wide_icc_data_contrasts$DrugvNeutral <- wide_icc_data_contrasts$drug - wide_icc_data_contrasts$neutral
  contrast_icc_data <- dcast(wide_icc_data_contrasts, id ~ time, value.var = 'DrugvNeutral')
  
  contrast_iccs <- ICC(contrast_icc_data[, 2:ncol(contrast_icc_data)])
  within_contrast_icc_intervals <- extract_icc_intervals(contrast_iccs)
  
  
  return(list(model = this_lme, plotframe = plot_frame, p = p, p2=p2, p3=p3, dset = long_data, within_visit_iccs = within_visit_iccs,
              within_iccs_simple_intervals = within_iccs_simple_intervals,
              within_contrast_icc_intervals = within_contrast_icc_intervals))
  
   
}


plot_one_beh <- function(measure, this_label, idps, limits_spag, limits_y = c(1, 4), prefix = 'dcr_', task = ''){
  #same as above, but for behavioral data
  #prefix will be dcr_tpre_ for tdcs/tacs and dcr_ for neurocaps 
  #measure will be box_rt, craving_rt, or craving
  #measure <- 'craving'
  #prefix <- 'dcr_tpre_'
  #this_label <- 'Craving'
  #idps <- idps_tdcs
  
  
  #idps_tacs$dcr_tpre_response_craving_[0-7]
  #idp_tdcs
  #dcr_tpre_response_box_rt_0
  #dcr_tpre_response_craving_rt_0
  #idps_neurocaps_ocr$dcr_response_box_craving[0-7]

  column_suffixes <- 0:7
  
  these_cols <- c('id', paste0(prefix, 'response_', measure, '_', column_suffixes))
  library(reshape2)
  one_dataset <- idps[, these_cols]
  long_data <- melt(one_dataset, id.vars = c('id'))
  long_data$variable <- as.character(long_data$variable)
  long_data$condition <- NA
  long_data$number <- substr(long_data$variable, nchar(long_data$variable), nchar(long_data$variable))
  long_data$condition[long_data$number %in% c('0', '2', '4', '6')] <- 'neutral'
  long_data$condition[long_data$number %in% c('1', '3', '5', '7')] <- 'drug'
  #put neutral in the intercept
  long_data$condition <- factor(long_data$condition, levels = c('neutral', 'drug'))
  #time = block number, just like for the imaging variables
  long_data$time <- NA
  long_data$time[long_data$number %in% c('0', '1')] <- 1
  long_data$time[long_data$number %in% c('2', '3')] <- 2
  long_data$time[long_data$number %in% c('4', '5')] <- 3
  long_data$time[long_data$number %in% c('6', '7')] <- 4

  #mean center on time
  long_data$time <- long_data$time - mean(long_data$time)

  library(lme4)
  library(lmerTest)
  library(ggplot2)
  library(sjstats) #for icc of mixed effects models
  #this_lme <- lmer(paste('value ~ condition * time + (1|id)'), data = long_data)
  #for checking for NA's, looks like it's all good now
  #print(long_data)
  this_lme <- lmer(paste('value ~ condition * time +  (1|id/condition)'), data = long_data)
  #print(summary(this_lme))
  within_visit_iccs <- sjstats::icc(this_lme)

  n_subjects <- length(unique(long_data$id))
  means <- aggregate(long_data$value, by = list(long_data$condition, long_data$time), FUN = mean, na.rm = TRUE)
  names(means)[1:2] <- c('condition', 'time')
  sds <- aggregate(long_data$value, by = list(long_data$condition, long_data$time), FUN = sd, na.rm = TRUE)
  names(sds)[1:2] <- c('condition', 'time')


  #print(this_label)
  mean_str <- paste(this_label, '_mean', sep = '')
  names(means)[3] <- mean_str 
  

  sd_str <- paste(this_label, '_sd', sep = '')
  se_str <- paste(this_label, '_se', sep = '')
  names(sds)[3] <- sd_str
  
  plot_frame <- merge(means, sds)
  plot_frame[, se_str] <- plot_frame[, sd_str] / sqrt(n_subjects)
  plot_frame$bar_top <- plot_frame[, mean_str] + plot_frame[, se_str]
  plot_frame$bar_bottom <- plot_frame[, mean_str] - plot_frame[, se_str]
  
  plot_frame$label <- this_label
  
  ebsize <- 1.5
  linesize <- 1
  #swap condition levels back for this plot, so that drug is red and neutral is blue
  plot_frame$condition <- factor(plot_frame$condition, levels = c('drug', 'neutral'))
  p <- ggplot(plot_frame) + geom_line(aes_string(x = 'time', y = mean_str, color = 'condition'), size = linesize) + geom_errorbar(aes_string(x = 'time', ymax = 'bar_top',                                                       ymin = 'bar_bottom', color = 'condition', width = 0.1), size = ebsize) +
                      theme(text = element_text(size=10)) + ggtitle(task) + ylab('') + xlab('') + ylim(limits_y)
  
  #print(p)
  
  
  p2 <- ggplot(data = long_data[long_data$condition == 'drug',], aes_string(x = 'time', y = 'value', group = 'id')) + 
    geom_line() + stat_summary(aes(group = 1), geom = "point", fun.y = mean, shape = 17, size = 3) + 
    stat_smooth(aes(group = 1)) + labs(
    x = "Time",
    y = this_label,
    title = "drug") 
  #print(p2)
  p3 <- ggplot(data = long_data[long_data$condition == 'neutral',], aes_string(x = 'time', y = 'value', group = 'id')) + 
    geom_line() + stat_summary(aes(group = 1), geom = "point", fun.y = mean, shape = 17, size = 3) + 
    stat_smooth(aes(group = 1)) + labs(
    x = "Time",
    y = this_label,
    title = "neutral") 
  #print(p3) 
  
  wide_icc_data <- dcast(long_data, id ~ variable)
  within_iccs_simple <- ICC(wide_icc_data[, 2:ncol(wide_icc_data)])
  within_iccs_simple_intervals <- extract_icc_intervals(within_iccs_simple)
  
  
  wide_icc_data_contrasts <- dcast(long_data, id + time ~ condition)
  wide_icc_data_contrasts$DrugvNeutral <- wide_icc_data_contrasts$drug - wide_icc_data_contrasts$neutral
  contrast_icc_data <- dcast(wide_icc_data_contrasts, id ~ time, value.var = 'DrugvNeutral')
  
  contrast_iccs <- ICC(contrast_icc_data[, 2:ncol(contrast_icc_data)])
  within_contrast_icc_intervals <- extract_icc_intervals(contrast_iccs)
  
  
  return(list(model = this_lme, plotframe = plot_frame, p = p, p2=p2, p3=p3, dset = long_data, within_visit_iccs = within_visit_iccs,
              within_iccs_simple_intervals = within_iccs_simple_intervals,
              within_contrast_icc_intervals = within_contrast_icc_intervals))
  
   
}


library(sjPlot)
## Learn more about sjPlot with 'browseVignettes("sjPlot")'.
library(grid)
library(gridExtra)
library(ggpubr)
## Loading required package: magrittr
library(psych)
## 
## Attaching package: 'psych'
## The following objects are masked from 'package:sjstats':
## 
##     pca, phi
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
#condition by time interaction
#1: VMPFC extracted at p < 0.001
#3: LSTG extracted at p < 0.001
#4: RSTG extracted at p < 0.001
#5: L Ventral Striatum extracted at p < 0.05 intersected with brainnetome ROIs 219 and 223
#6: R Ventral Striatum extracted at p < 0.05 intersected with brainnetome ROIs 220 and 224
#8: R Amygdala extracted at p < 0.05 intersected with brainnetome ROIs 212 and 214

#main effect of condition
#9: LIFG extracted at p < 0.001
#10: RIFG at p < 0.001
#12: RDLPFC at p < 0.001

extract_icc_intervals <- function(this_icc){
  
  single_absolute <- round(this_icc$results$ICC[1], 3)
  single_absolute_lower <- this_icc$results$`lower bound`[1]
  single_absolute_upper <- this_icc$results$`upper bound`[1]
  single_fixed <- round(this_icc$results$ICC[3], 3)
  single_fixed_lower <- this_icc$results$`lower bound`[3]
  single_fixed_upper <- this_icc$results$`upper bound`[3]
  average_fixed <- round(this_icc$results$ICC[6], 3)
  average_fixed_lower <- this_icc$results$`lower bound`[6]
  average_fixed_upper <- this_icc$results$`upper bound`[6]
  this_row <- data.frame(t(c(ICC1 = single_absolute, ICC1_lower = single_absolute_lower, ICC1_upper = single_absolute_upper, ICC3 = single_fixed,
                             ICC3_lower = single_fixed_lower, ICC3_upper = single_fixed_upper,
                             ICC3k = average_fixed, ICC3k_lower = average_fixed_lower, ICC3k_upper = average_fixed_upper)))
  return(this_row)   
}

compare_models <- function(tdcs_model, neurocaps_mcr_model, neurocaps_ocr_model, tacs_model){
  #prints off tables of z-scores and p-values for the coefficients of condition and condition*time compared between all pairs of models
  #tdcs_summary <- summary(tdcs_model)
  #fill this with z-scores testing for differences between coefficients of different models
  ztable_condition <- data.frame(tdcs=c(NA,NA,NA,NA), neurocaps_mcr = c(NA,NA,NA,NA),
                       neurocaps_ocr = c(NA,NA,NA,NA), tacs=c(NA,NA,NA,NA))
  rownames(ztable_condition) <- c('tdcs', 'neurocaps_mcr', 'neurocaps_ocr', 'tacs')
  
  ztable_conditionbytime <- data.frame(tdcs=c(NA,NA,NA,NA), neurocaps_mcr = c(NA,NA,NA,NA),
                       neurocaps_ocr = c(NA,NA,NA,NA), tacs=c(NA,NA,NA,NA))
  rownames(ztable_conditionbytime) <- c('tdcs', 'neurocaps_mcr', 'neurocaps_ocr', 'tacs')
  
  all_models <- list(summary(tdcs_model), summary(neurocaps_mcr_model), summary(neurocaps_ocr_model), summary(tacs_model))
  for (i in 1:3){
    for (j in (i+1):4){
      #condition beta
      beta1 <- all_models[[i]]$coefficients[2,1]
      #condition SE
      se1 <- all_models[[i]]$coefficients[2,2]
      #condition beta
      beta2 <- all_models[[j]]$coefficients[2,1]
      #condition SE
      se2 <- all_models[[j]]$coefficients[2,2]
      this_z <- (beta1 - beta2) / sqrt(se1 * se1 + se2 * se2)
      ztable_condition[i,j] <- ztable_condition[j,i] <- this_z
      
      
      #conditionbytime beta--row 5 for imaging and 4 for behavioral (since motion is not included there)--but always the last row
      amodel <- all_models[[i]]$coefficients
      cbytime_idx <- nrow(amodel)
      beta1 <- all_models[[i]]$coefficients[cbytime_idx,1]
      #conditionbytime SE
      se1 <- all_models[[i]]$coefficients[cbytime_idx,2]
      #conditionbytime beta
      beta2 <- all_models[[j]]$coefficients[cbytime_idx,1]
      #conditionbytime SE
      se2 <- all_models[[j]]$coefficients[cbytime_idx,2]
      this_z <- (beta1 - beta2) / sqrt(se1 * se1 + se2 * se2)
      ztable_conditionbytime[i,j] <- ztable_conditionbytime[j,i] <- this_z
      
    }
  }
  
  ptable_condition <- ztable_condition
  ptable_conditionbytime <- ztable_conditionbytime
  for (n in names(ptable_condition)){
    ptable_condition[,n] <- 2*pnorm(-abs(ptable_condition[,n]))
    ptable_conditionbytime[,n] <- 2*pnorm(-abs(ptable_conditionbytime[,n]))
  }
  print('###Condition z-scores###')
  print(ztable_condition)
  print('###Condition by time z-scores###')
  print(ztable_conditionbytime)
    
  print('###Condition p-values###')
  print(ptable_condition)
  print('###Condition by time p-values###')
  print(ptable_conditionbytime)
  
  
  
  
  
}
  


plot_combination <- function(roi, label, limits_spag = c(-2,2), limits_y = c(-1, 1), imaging = TRUE, forest_range = c(-0.5, 0.5)){

  #roi <- '9'
  #label <- 'LIFG'
  #limits_y = c(-0.3, 0.5)
  #limits_spag = c(-2, 2)
  
  if (imaging){
    tdcs_list <- plot_one(roi, label, idps_tdcs, limits_spag = limits_spag, limits_y = limits_y, prefix = 'dcr_tpre_', task = '1212MCR')
    neurocaps_mcr_list <- plot_one(roi, label, idps_neurocaps_mcr, limits_spag = limits_spag, limits_y = limits_y, prefix = 'dcr_', task = 'WiRMCR')
    neurocaps_ocr_list <- plot_one(roi, label, idps_neurocaps_ocr, limits_spag = limits_spag, limits_y = limits_y, prefix = 'dcr_', task = '1212OCR')
    tacs_list <- plot_one(roi, label, idps_tacs, limits_spag = limits_spag, limits_y = limits_y, prefix = 'dcr_tpre_', task = '1212OCR2')
  } else {
    tdcs_list <- plot_one_beh(roi, label, idps_tdcs, limits_spag = limits_spag, limits_y = limits_y, prefix = 'dcr_tpre_', task = '1212MCR')
    neurocaps_mcr_list <- plot_one_beh(roi, label, idps_neurocaps_mcr, limits_spag = limits_spag, limits_y = limits_y, prefix = 'dcr_', task = 'WiRMCR')
    neurocaps_ocr_list <- plot_one_beh(roi, label, idps_neurocaps_ocr, limits_spag = limits_spag, limits_y = limits_y, prefix = 'dcr_', task = '1212OCR')
    tacs_list <- plot_one_beh(roi, label, idps_tacs, limits_spag = limits_spag, limits_y = limits_y, prefix = 'dcr_tpre_', task = '1212OCR2')
  }
  compare_models(tdcs_list$model, neurocaps_mcr_list$model, neurocaps_ocr_list$model, tacs_list$model)
  
  
  #print model with restricted range
  print(plot_models(tdcs_list$model, neurocaps_mcr_list$model, neurocaps_ocr_list$model, tacs_list$model, m.labels = c('tDCS', 'WirMCR', '1212OCR', '1212OCR2')) + ylim(forest_range))
  #plot model without restricting range so you can see where 'motion' falls
  print(plot_models(tdcs_list$model, neurocaps_mcr_list$model, neurocaps_ocr_list$model, tacs_list$model, m.labels = c('tDCS', 'WirMCR', '1212OCR', '1212OCR2')))# + ylim(-0.5, 0.5))
  
  #from:
  #https://stats.stackexchange.com/questions/93540/testing-equality-of-coefficients-from-two-different-regressions
  #to test for a difference in betas, compute:
  #Z = (beta1 - beta2)/sqrt(SE(beta1)^2 + SE(beta2)^2)
  pwidth = 6
  pheight = 8
  print(ggarrange(tdcs_list$p, neurocaps_mcr_list$p, neurocaps_ocr_list$p, tacs_list$p, ncol = 2, nrow = 2, common.legend = TRUE, legend = 'bottom'))
  ggsave(paste0(label, '-lineplot.png') , ggarrange(tdcs_list$p, neurocaps_mcr_list$p, neurocaps_ocr_list$p, tacs_list$p, ncol = 2, nrow = 2, 
                                                    common.legend = TRUE, legend = 'bottom'),
         width = pwidth, height = pheight)
  print(ggarrange(tdcs_list$p2, neurocaps_mcr_list$p2, neurocaps_ocr_list$p2, tacs_list$p2, ncol = 2, nrow = 2, common.legend = TRUE, legend = 'bottom'))
  
  print(ggarrange(tdcs_list$p3, neurocaps_mcr_list$p3, neurocaps_ocr_list$p3, tacs_list$p3, ncol = 2, nrow = 2, common.legend = TRUE, legend = 'bottom'))
  
  
  #get ICCs for each of the 8 conditions
  icc_data <- merge(neurocaps_ocr_list$dset, tacs_list$dset, by = c('id', 'time', 'condition'))
  all_iccs <- NULL
  mean_values <- NA
  for (con in unique(icc_data$condition)){
    for (tim in unique(icc_data$time)){
      
      for_icc <- icc_data[icc_data$time == tim & icc_data$condition == con, c('value.x', 'value.y')]
      this_icc <- ICC(for_icc)
      this_row <- extract_icc_intervals(this_icc)
      this_row$var <- paste(tim, con)      
      all_iccs <- rbind(all_iccs, this_row)
      
      #get values for drug and neutral at this time point to put into mean_values frame
      for_means <- icc_data[icc_data$time == tim & icc_data$condition == con, c('id', 'value.x', 'value.y')]
      value_cols <- paste0(con, tim, c('NeuroCaps', 'tACS'))
      names(for_means) <- c('id', value_cols)
      mean_values <- merge(mean_values, for_means, all = TRUE)
      
    }
    #get ICCs for mean of drug and neutral
    for_icc <- icc_data[icc_data$condition == con, c('id', 'value.x', 'value.y')]
    for_icc <- aggregate(for_icc[, c('value.x', 'value.y')], by = list(for_icc$id), FUN = mean)
    value_cols <- paste0(con, c('NeuroCaps', 'tACS'))
    names(for_icc) <- c('id', value_cols)
    mean_values <- merge(mean_values, for_icc, all = TRUE)
    this_icc <- ICC(for_icc[, value_cols])
    this_row <- extract_icc_intervals(this_icc)
    this_row$var <- paste(con)      
    all_iccs <- rbind(all_iccs, this_row)
  }
  
  for (tim in unique(icc_data$time)){
    neurocaps_colname <- paste0('NeuroCapsDvN', tim)
    tacs_colname <- paste0('tACSDvN', tim)
    mean_values[, neurocaps_colname] <- mean_values[, paste0('drug', tim, 'NeuroCaps')] - mean_values[, paste0('neutral', tim, 'NeuroCaps')]
    mean_values[, tacs_colname] <- mean_values[, paste0('drug', tim, 'tACS')] - mean_values[, paste0('neutral', tim, 'tACS')]
    this_icc <- ICC(mean_values[, c(neurocaps_colname, tacs_colname)])
    this_row <- extract_icc_intervals(this_icc)
    this_row$var <- paste0(tim, 'DrugvNeutral')      
    all_iccs <- rbind(all_iccs, this_row)
  }

  #get early and late mean values
  mean_values$NeuroCapsDvNEarly <- (mean_values$`NeuroCapsDvN-1.5` + mean_values$`NeuroCapsDvN-0.5`) / 2
  mean_values$NeuroCapsDvNLate <- (mean_values$`NeuroCapsDvN1.5` + mean_values$`NeuroCapsDvN0.5`) / 2
  
  mean_values$tACSDvNEarly <- (mean_values$`tACSDvN-1.5` + mean_values$`tACSDvN-0.5`) / 2
  mean_values$tACSDvNLate <- (mean_values$`tACSDvN1.5` + mean_values$`tACSDvN0.5`) / 2
  
  #compute early and late ICCs
  this_icc <- ICC(mean_values[, c('NeuroCapsDvNEarly', 'tACSDvNEarly')])
  this_row <- extract_icc_intervals(this_icc)
  this_row$var <- paste('DrugvNeutralEarly')      
  all_iccs <- rbind(all_iccs, this_row)
  
  this_icc <- ICC(mean_values[, c('NeuroCapsDvNLate', 'tACSDvNLate')])
  this_row <- extract_icc_intervals(this_icc)
  this_row$var <- paste('DrugvNeutralLate')      
  all_iccs <- rbind(all_iccs, this_row)
  
  
  
  mean_values$NeuroCapsDvN <- mean_values$drugNeuroCaps - mean_values$neutralNeuroCaps
  mean_values$tACSDvN <- mean_values$drugtACS - mean_values$neutraltACS
  
  this_icc <- ICC(mean_values[, c('NeuroCapsDvN', 'tACSDvN')])
  this_row <- extract_icc_intervals(this_icc)
  this_row$var <- paste('DrugvNeutral')      
  all_iccs <- rbind(all_iccs, this_row)
  
  #get ICCs for mean contrast of drug and neutral in each block and in blocks 1/2 and 3/4 separately
  
  
  
  melted_iccs <- all_iccs[, c('var', 'ICC1', 'ICC1_lower', 'ICC1_upper')]
  names(melted_iccs) <- c('var', 'ICC', 'lower', 'upper')
  melted_iccs$variable <- 'ICC1'

  melted_iccs3 <- all_iccs[, c('var', 'ICC3', 'ICC3_lower', 'ICC3_upper')]
  names(melted_iccs3) <- c('var', 'ICC', 'lower', 'upper')
  melted_iccs3$variable <- 'ICC3'
  
  melted_iccs <- rbind(melted_iccs, melted_iccs3)
  #melted_iccs$var <- factor(melted_iccs$var, levels = rev(to_plot))
  melted_iccs$var <- factor(melted_iccs$var, levels = c('DrugvNeutral', 'drug', 'neutral',
                                                        "DrugvNeutralLate", "DrugvNeutralEarly",
                                                        "-1.5DrugvNeutral", "-0.5DrugvNeutral",  "0.5DrugvNeutral",   "1.5DrugvNeutral",
                                                           '1.5 drug', '0.5 drug', '-0.5 drug', '-1.5 drug',
                                                           '1.5 neutral', '0.5 neutral', '-0.5 neutral', '-1.5 neutral'))
  p <- ggplot(melted_iccs) + geom_bar(aes(x = var, fill = variable, y = ICC), stat = 'identity', position = 'dodge') + coord_flip() +
    xlab('Variable') + ylab('ICC') + ggtitle(paste0('Between Session: ', label)) + geom_errorbar(aes(x = var, fill = variable, ymin = lower, ymax = upper), stat = 'identity', position = 'dodge')
  print(p)
  print(all_iccs)
  ggsave(paste0('ICCs_', label, '.png'), p)
  
  
  #plot ICCs that come from the LME
  within_icc_results <- as.data.frame(rbind(tdcs_list$within_visit_iccs, neurocaps_mcr_list$within_visit_iccs, neurocaps_ocr_list$within_visit_iccs, tacs_list$within_visit_iccs))
  within_icc_results$dset <- c('tDCS', 'WirMCR', '1212OCR', '1212OCR2')
  within_icc_results$id_condition <- within_icc_results$id + within_icc_results$`condition:id`
  melted_within_data <- melt(within_icc_results, id.vars = c('dset', 'condition:id'))
  melted_within_data$dset <- factor(melted_within_data$dset, levels = c('1212OCR',  '1212OCR2','WirMCR', 'tDCS'))
  p <- ggplot(data = melted_within_data) + geom_bar(aes(x = dset, y = value, color = variable, fill = variable), stat = 'identity', position = 'dodge') +
    coord_flip() + ggtitle(paste0('Within Session ICCs from LME for: ', label))
  print(p)
  
  #plot within-session ICCs from psych's ICC function, treating all conditions as the same (should be similar to 'id' above)
  tdcs_list$within_iccs_simple_intervals$var <- 'tDCS'
  all_iccs <- tdcs_list$within_iccs_simple_intervals
  neurocaps_mcr_list$within_iccs_simple_intervals$var <- 'WirMCR'
  all_iccs <- rbind(all_iccs, neurocaps_mcr_list$within_iccs_simple_intervals)
  neurocaps_ocr_list$within_iccs_simple_intervals$var <- '1212OCR'
  all_iccs <- rbind(all_iccs, neurocaps_ocr_list$within_iccs_simple_intervals)
  tacs_list$within_iccs_simple_intervals$var <- '1212OCR2'
  all_iccs <- rbind(all_iccs, tacs_list$within_iccs_simple_intervals)
  
  ###plot--should functionalize this
  melted_iccs <- all_iccs[, c('var', 'ICC1', 'ICC1_lower', 'ICC1_upper')]
  names(melted_iccs) <- c('var', 'ICC', 'lower', 'upper')
  melted_iccs$variable <- 'ICC1'

  melted_iccs3 <- all_iccs[, c('var', 'ICC3', 'ICC3_lower', 'ICC3_upper')]
  names(melted_iccs3) <- c('var', 'ICC', 'lower', 'upper')
  melted_iccs3$variable <- 'ICC3'
  
  melted_iccs3k <- all_iccs[, c('var', 'ICC3k', 'ICC3k_lower', 'ICC3k_upper')]
  names(melted_iccs3k) <- c('var', 'ICC', 'lower', 'upper')
  melted_iccs3k$variable <- 'ICC3k'
  
  melted_iccs <- rbind(melted_iccs, melted_iccs3)
  melted_iccs <- rbind(melted_iccs, melted_iccs3k)
  #melted_iccs$var <- factor(melted_iccs$var, levels = rev(to_plot))
  melted_iccs$var <- factor(melted_iccs$var, levels = c('1212OCR2',  '1212OCR','WirMCR', 'tDCS'))
  p <- ggplot(melted_iccs) + geom_bar(aes(x = var, fill = variable, y = ICC), stat = 'identity', position = 'dodge') + coord_flip() +
    xlab('Variable') + ylab('ICC') + ggtitle(paste0('Within Session Simple ICC: ', label)) + geom_errorbar(aes(x = var, fill = variable, ymin = lower, ymax = upper), stat = 'identity', position = 'dodge')
  print(p)
  
  
  #plot within-session ICCs from psych's ICC function, using contrasts
  tdcs_list$within_contrast_icc_intervals$var <- 'tDCS'
  all_iccs <- tdcs_list$within_contrast_icc_intervals
  neurocaps_mcr_list$within_contrast_icc_intervals$var <- 'WirMCR'
  all_iccs <- rbind(all_iccs, neurocaps_mcr_list$within_contrast_icc_intervals)
  neurocaps_ocr_list$within_contrast_icc_intervals$var <- '1212OCR'
  all_iccs <- rbind(all_iccs, neurocaps_ocr_list$within_contrast_icc_intervals)
  tacs_list$within_contrast_icc_intervals$var <- '1212OCR2'
  all_iccs <- rbind(all_iccs, tacs_list$within_contrast_icc_intervals)
  
  ###plot--should functionalize this
  melted_iccs <- all_iccs[, c('var', 'ICC1', 'ICC1_lower', 'ICC1_upper')]
  names(melted_iccs) <- c('var', 'ICC', 'lower', 'upper')
  melted_iccs$variable <- 'ICC1'

  melted_iccs3 <- all_iccs[, c('var', 'ICC3', 'ICC3_lower', 'ICC3_upper')]
  names(melted_iccs3) <- c('var', 'ICC', 'lower', 'upper')
  melted_iccs3$variable <- 'ICC3'
  
  melted_iccs3k <- all_iccs[, c('var', 'ICC3k', 'ICC3k_lower', 'ICC3k_upper')]
  names(melted_iccs3k) <- c('var', 'ICC', 'lower', 'upper')
  melted_iccs3k$variable <- 'ICC3k'
  
  melted_iccs <- rbind(melted_iccs, melted_iccs3)
  melted_iccs <- rbind(melted_iccs, melted_iccs3k)
  #melted_iccs$var <- factor(melted_iccs$var, levels = rev(to_plot))
  melted_iccs$var <- factor(melted_iccs$var, levels = c('1212OCR2',  '1212OCR','WirMCR', 'tDCS'))
  p <- ggplot(melted_iccs) + geom_bar(aes(x = var, fill = variable, y = ICC), stat = 'identity', position = 'dodge') + coord_flip() +
    xlab('Variable') + ylab('ICC') + ggtitle(paste0('Within Session Contrast ICC: ', label)) + geom_errorbar(aes(x = var, fill = variable, ymin = lower, ymax = upper), stat = 'identity', position = 'dodge')
  print(p)
  
  ###
  
  
  
  
  print('###tDCS MCR###')
  print(summary(tdcs_list$model))
  #print(anova(tdcs_list$model, type = 'marginal'))
  print('###Neurocaps MCR###')
  print(summary(neurocaps_mcr_list$model))
  print('###Neurocaps OCR###')
  print(summary(neurocaps_ocr_list$model))
  print('###tACS OCR###')
  print(summary(tacs_list$model))
  
  
  ##make model output table to save
  tablefile <- paste0(label, '-modeltable.csv')
  write.table('tDCS Model Output', tablefile)
  write.table(round(summary(tdcs_list$model)$coefficients, digits = 3), file = tablefile,
            row.names = TRUE, append = TRUE, sep = ',')
  
  write.table('Neurocaps MCR Model Output', tablefile, append = TRUE)
  write.table(round(summary(neurocaps_mcr_list$model)$coefficients, digits = 3), file = tablefile,
            row.names = TRUE, append = TRUE, sep = ',')
  
  write.table('Neurocaps OCR Model Output', tablefile, append = TRUE)
  write.table(round(summary(neurocaps_ocr_list$model)$coefficients, digits = 3), file = tablefile,
            row.names = TRUE, append = TRUE, sep = ',')
  
  write.table('tACS OCR Model Output', tablefile, append = TRUE)
  write.table(round(summary(tacs_list$model)$coefficients, digits = 3), file = tablefile,
            row.names = TRUE, append = TRUE, sep = ',')
  
  
}

Condition

LIFG

knitr::include_graphics("LIFG_z16.png")

plot_combination('9', 'LIFG', limits_y = c(-0.15, 0.45))
## [1] "###Condition z-scores###"
##                    tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                 NA     1.1616086     0.8054013  0.3833178
## neurocaps_mcr 1.1616086            NA    -0.1122026 -0.5502508
## neurocaps_ocr 0.8054013    -0.1122026            NA -0.3762163
## tacs          0.3833178    -0.5502508    -0.3762163         NA
## [1] "###Condition by time z-scores###"
##                     tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                  NA      1.222763    -0.3421076 -1.831715
## neurocaps_mcr  1.2227631            NA    -1.2144775 -2.462810
## neurocaps_ocr -0.3421076     -1.214477            NA -1.190708
## tacs          -1.8317149     -2.462810    -1.1907079        NA
## [1] "###Condition p-values###"
##                    tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                 NA     0.2453945     0.4205881 0.7014841
## neurocaps_mcr 0.2453945            NA     0.9106628 0.5821474
## neurocaps_ocr 0.4205881     0.9106628            NA 0.7067561
## tacs          0.7014841     0.5821474     0.7067561        NA
## [1] "###Condition by time p-values###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA    0.22141917     0.7322699 0.06699391
## neurocaps_mcr 0.22141917            NA     0.2245654 0.01378528
## neurocaps_ocr 0.73226988    0.22456544            NA 0.23376828
## tacs          0.06699391    0.01378528     0.2337683         NA
## Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

##      ICC1  ICC1_lower ICC1_upper  ICC3  ICC3_lower ICC3_upper ICC3k ICC3k_lower ICC3k_upper               var
## 1  -0.011 -0.41613414  0.4015172 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192         -0.5 drug
## 2   0.570  0.21179397  0.7944981 0.570  0.20462615  0.7956329 0.726  0.33973387   0.8861866         -1.5 drug
## 3   0.372 -0.04089754  0.6791548 0.385 -0.03340966  0.6887739 0.556 -0.06912889   0.8157088          0.5 drug
## 4   0.203 -0.22190266  0.5668081 0.203 -0.22901444  0.5688981 0.338 -0.59408230   0.7252200          1.5 drug
## 5   0.535  0.16369113  0.7753645 0.543  0.16742384  0.7810495 0.704  0.28682614   0.8770666              drug
## 6   0.060 -0.35530845  0.4597212 0.089 -0.33631783  0.4846977 0.164 -1.01349062   0.6529244      -0.5 neutral
## 7   0.383 -0.02882437  0.6856138 0.387 -0.03072059  0.6901861 0.559 -0.06338850   0.8166983      -1.5 neutral
## 8   0.069 -0.34750317  0.4667153 0.090 -0.33572472  0.4852091 0.165 -1.01079998   0.6533882       0.5 neutral
## 9   0.031 -0.38099588  0.4359664 0.066 -0.35686693  0.4666275 0.124 -1.10977631   0.6363272       1.5 neutral
## 10  0.219 -0.20686902  0.5774081 0.219 -0.21402982  0.5794607 0.359 -0.54462579   0.7337450           neutral
## 11  0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192  -0.5DrugvNeutral
## 12  0.243 -0.18220152  0.5942463 0.243 -0.18943545  0.5962381 0.391 -0.46741607   0.7470541  -1.5DrugvNeutral
## 13  0.511  0.13178270  0.7620211 0.511  0.12441301  0.7633116 0.677  0.22129414   0.8657705   0.5DrugvNeutral
## 14 -0.032 -0.43372330  0.3833580 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192   1.5DrugvNeutral
## 15  0.238 -0.18665261  0.5912574 0.238 -0.19387401  0.5932602 0.385 -0.48100177   0.7447122 DrugvNeutralEarly
## 16  0.327 -0.09261566  0.6501636 0.329 -0.09801728  0.6531131 0.495 -0.21733738   0.7901614  DrugvNeutralLate
## 17  0.491  0.10440168  0.7501278 0.491  0.09698501  0.7514737 0.658  0.17682102   0.8581045      DrugvNeutral
## Saving 7 x 5 in image

## [1] "###tDCS MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 125.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9228 -0.5555  0.0131  0.5274  3.4993 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.00000  0.0000  
##  id           (Intercept) 0.01771  0.1331  
##  Residual                 0.06172  0.2484  
## Number of obs: 520, groups:  condition:id, 130; id, 65
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)          0.232357   0.042932  71.923445   5.412 7.80e-07 ***
## conditiondrug        0.140665   0.021789 452.000000   6.456 2.78e-10 ***
## time                 0.013784   0.013781 452.000000   1.000  0.31774    
## motion              -1.119903   0.356675  63.000000  -3.140  0.00257 ** 
## conditiondrug:time  -0.003891   0.019489 452.000000  -0.200  0.84186    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.254                     
## time         0.000  0.000              
## motion      -0.851  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000
## [1] "###Neurocaps MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 142.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.8048 -0.4694 -0.0166  0.5285  3.8206 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.005082 0.07129 
##  id           (Intercept) 0.022843 0.15114 
##  Residual                 0.084911 0.29139 
## Number of obs: 232, groups:  condition:id, 58; id, 29
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)  
## (Intercept)          0.05421    0.08188  30.90438   0.662   0.5128  
## conditiondrug        0.08509    0.04260  28.00000   1.998   0.0556 .
## time                 0.02114    0.02420 172.00000   0.874   0.3836  
## motion              -0.13958    0.71730  27.00000  -0.195   0.8472  
## conditiondrug:time  -0.05205    0.03422 172.00000  -1.521   0.1301  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.260                     
## time         0.000  0.000              
## motion      -0.864  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000
## [1] "###Neurocaps OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 106.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7934 -0.5240  0.0839  0.5783  2.5766 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.012444 0.11155 
##  id           (Intercept) 0.007477 0.08647 
##  Residual                 0.084054 0.28992 
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)
## (Intercept)          0.07747    0.07432  26.25739   1.042    0.307
## conditiondrug        0.09291    0.05515  21.00000   1.685    0.107
## time                 0.02024    0.02764 130.00000   0.732    0.465
## motion               0.30032    0.56748  20.00000   0.529    0.602
## conditiondrug:time   0.01105    0.03909 130.00000   0.283    0.778
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.371                     
## time         0.000  0.000              
## motion      -0.814  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000
## [1] "###tACS OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 116.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.1737 -0.5107 -0.0123  0.4421  2.9365 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.003373 0.05808 
##  id           (Intercept) 0.028954 0.17016 
##  Residual                 0.088214 0.29701 
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)  
## (Intercept)          0.01767    0.10289  22.30565   0.172   0.8652  
## conditiondrug        0.12043    0.04808  21.00000   2.505   0.0206 *
## time                -0.05792    0.02832 130.00000  -2.045   0.0428 *
## motion               1.37727    1.00953  20.00000   1.364   0.1876  
## conditiondrug:time   0.07769    0.04005 130.00000   1.940   0.0546 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.234                     
## time         0.000  0.000              
## motion      -0.875  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000

RIFG

knitr::include_graphics("RIFG_z7.png")

plot_combination('10', 'RIFG', limits_y = c(-0.15, 0.5))
## [1] "###Condition z-scores###"
##                    tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                 NA     0.7123509     0.7960760 1.4727623
## neurocaps_mcr 0.7123509            NA     0.1962350 0.7743601
## neurocaps_ocr 0.7960760     0.1962350            NA 0.5133078
## tacs          1.4727623     0.7743601     0.5133078        NA
## [1] "###Condition by time z-scores###"
##                     tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                  NA    0.86275495    0.95744543 -1.402705
## neurocaps_mcr  0.8627550            NA    0.09483428 -1.808677
## neurocaps_ocr  0.9574454    0.09483428            NA -1.872828
## tacs          -1.4027053   -1.80867704   -1.87282849        NA
## [1] "###Condition p-values###"
##                    tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                 NA     0.4762475     0.4259879 0.1408151
## neurocaps_mcr 0.4762475            NA     0.8444263 0.4387179
## neurocaps_ocr 0.4259879     0.8444263            NA 0.6077361
## tacs          0.1408151     0.4387179     0.6077361        NA
## [1] "###Condition by time p-values###"
##                    tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                 NA    0.38827221    0.33834248 0.16070473
## neurocaps_mcr 0.3882722            NA    0.92444646 0.07050119
## neurocaps_ocr 0.3383425    0.92444646            NA 0.06109207
## tacs          0.1607047    0.07050119    0.06109207         NA
## Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

##      ICC1  ICC1_lower ICC1_upper  ICC3  ICC3_lower ICC3_upper ICC3k ICC3k_lower ICC3k_upper               var
## 1   0.133 -0.28977837  0.5154618 0.133 -0.29662688  0.5177233 0.235 -0.84344104   0.6822367         -0.5 drug
## 2   0.461  0.06586948  0.7326521 0.550  0.17634251  0.7846078 0.709  0.29981490   0.8793056         -1.5 drug
## 3   0.600  0.25478919  0.8106739 0.600  0.24776965  0.8117287 0.750  0.39714005   0.8960820          0.5 drug
## 4   0.140 -0.28298280  0.5208768 0.140 -0.28986075  0.5231210 0.246 -0.81634907   0.6869067          1.5 drug
## 5   0.472  0.08058437  0.7394312 0.472  0.07313599  0.7408262 0.641  0.13630331   0.8511202              drug
## 6   0.282 -0.14080887  0.6210504 0.282 -0.14814511  0.6229417 0.440 -0.34781773   0.7676698      -0.5 neutral
## 7   0.358 -0.05764855  0.6700045 0.358 -0.06511296  0.6717010 0.527 -0.13929590   0.8036138      -1.5 neutral
## 8   0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192       0.5 neutral
## 9  -0.045 -0.44361712  0.3728581 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192       1.5 neutral
## 10  0.063 -0.35274973  0.4620252 0.063 -0.35929266  0.4644480 0.119 -1.12155001   0.6342977           neutral
## 11  0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192  -0.5DrugvNeutral
## 12  0.249 -0.17594668  0.5984104 0.374 -0.04657415  0.6817805 0.544 -0.09769852   0.8107842  -1.5DrugvNeutral
## 13  0.608  0.26729981  0.8152265 0.608  0.26032865  0.8162586 0.756  0.41311232   0.8988352   0.5DrugvNeutral
## 14 -0.011 -0.41595961  0.4016943 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192   1.5DrugvNeutral
## 15  0.112 -0.30899658  0.4997915 0.134 -0.29592416  0.5182869 0.236 -0.84060310   0.6827259 DrugvNeutralEarly
## 16  0.650  0.33015210  0.8371177 0.650  0.32345970  0.8380384 0.788  0.48880929   0.9118835  DrugvNeutralLate
## 17  0.499  0.11567526  0.7550758 0.499  0.10827656  0.7563987 0.666  0.19539629   0.8613064      DrugvNeutral
## Saving 7 x 5 in image

## [1] "###tDCS MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 127.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2004 -0.5378  0.0284  0.5999  3.5567 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.00000  0.0000  
##  id           (Intercept) 0.01229  0.1108  
##  Residual                 0.06403  0.2530  
## Number of obs: 520, groups:  condition:id, 130; id, 65
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)         1.272e-01  3.872e-02  7.469e+01   3.284  0.00156 ** 
## conditiondrug       1.661e-01  2.219e-02  4.520e+02   7.485 3.76e-13 ***
## time                4.917e-03  1.404e-02  4.520e+02   0.350  0.72625    
## motion             -5.992e-01  3.186e-01  6.300e+01  -1.880  0.06468 .  
## conditiondrug:time -6.029e-04  1.985e-02  4.520e+02  -0.030  0.97578    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.287                     
## time         0.000  0.000              
## motion      -0.842  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000
## [1] "###Neurocaps MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 179.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.0111 -0.3894  0.0024  0.4082  5.8371 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.00122  0.03493 
##  id           (Intercept) 0.02297  0.15155 
##  Residual                 0.10432  0.32299 
## Number of obs: 232, groups:  condition:id, 58; id, 29
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)          0.18159    0.08263  30.98218   2.198  0.03558 * 
## conditiondrug        0.13141    0.04339  28.00000   3.029  0.00523 **
## time                 0.01599    0.02682 172.00001   0.596  0.55194   
## motion              -0.35930    0.72344  27.00000  -0.497  0.62345   
## conditiondrug:time  -0.03754    0.03793 172.00001  -0.990  0.32374   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.263                     
## time         0.000  0.000              
## motion      -0.864  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000
## [1] "###Neurocaps OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 107
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4826 -0.5020 -0.0037  0.5033  4.3124 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.014729 0.12136 
##  id           (Intercept) 0.004287 0.06548 
##  Residual                 0.084736 0.29110 
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)  
## (Intercept)          0.09665    0.07207  27.25103   1.341   0.1909  
## conditiondrug        0.11733    0.05714  21.00000   2.053   0.0527 .
## time                 0.02199    0.02775 130.00000   0.792   0.4297  
## motion               0.11276    0.54399  20.00000   0.207   0.8379  
## conditiondrug:time  -0.04272    0.03925 130.00000  -1.088   0.2785  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.396                     
## time         0.000  0.000              
## motion      -0.805  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000
## [1] "###tACS OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 111.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5691 -0.4648 -0.0028  0.4447  4.3509 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.01733  0.1317  
##  id           (Intercept) 0.04992  0.2234  
##  Residual                 0.07431  0.2726  
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)
## (Intercept)         -0.02137    0.13084  22.00018  -0.163    0.872
## conditiondrug        0.07586    0.05714  21.00000   1.328    0.199
## time                -0.03981    0.02599 130.00000  -1.532    0.128
## motion               1.22228    1.28850  20.00000   0.949    0.354
## conditiondrug:time   0.05800    0.03676 130.00000   1.578    0.117
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.218                     
## time         0.000  0.000              
## motion      -0.879  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000

RDLPFC

knitr::include_graphics("RDLPFC_z7.png")

plot_combination('12', 'RDLPFC', limits_y = c(-0.45, 0.1))
## [1] "###Condition z-scores###"
##                      tdcs neurocaps_mcr neurocaps_ocr        tacs
## tdcs                   NA    0.31425131   0.009064596  0.21861556
## neurocaps_mcr 0.314251312            NA  -0.212396545 -0.05095833
## neurocaps_ocr 0.009064596   -0.21239655            NA  0.15452834
## tacs          0.218615559   -0.05095833   0.154528341          NA
## [1] "###Condition by time z-scores###"
##                      tdcs neurocaps_mcr neurocaps_ocr        tacs
## tdcs                   NA    0.05074138    0.02546767 -0.07062211
## neurocaps_mcr  0.05074138            NA   -0.02233374 -0.09538805
## neurocaps_ocr  0.02546767   -0.02233374            NA -0.07761132
## tacs          -0.07062211   -0.09538805   -0.07761132          NA
## [1] "###Condition p-values###"
##                    tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                 NA     0.7533302     0.9927676 0.8269495
## neurocaps_mcr 0.7533302            NA     0.8317977 0.9593587
## neurocaps_ocr 0.9927676     0.8317977            NA 0.8771932
## tacs          0.8269495     0.9593587     0.8771932        NA
## [1] "###Condition by time p-values###"
##                    tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                 NA     0.9595316     0.9796819 0.9436985
## neurocaps_mcr 0.9595316            NA     0.9821817 0.9240066
## neurocaps_ocr 0.9796819     0.9821817            NA 0.9381372
## tacs          0.9436985     0.9240066     0.9381372        NA
## Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

##     ICC1  ICC1_lower ICC1_upper  ICC3  ICC3_lower ICC3_upper ICC3k ICC3k_lower ICC3k_upper               var
## 1  0.037 -0.37551627  0.4411311 0.037 -0.38193421  0.4436122 0.072 -1.23590147   0.6145864         -0.5 drug
## 2  0.153 -0.27069876  0.5305039 0.202 -0.23015741  0.5680815 0.337 -0.59793370   0.7245561         -1.5 drug
## 3  0.071 -0.34593490  0.4681084 0.071 -0.35251375  0.4705138 0.133 -1.08886868   0.6399311          0.5 drug
## 4  0.287 -0.13637989  0.6238166 0.300 -0.12966353  0.6343391 0.461 -0.29796184   0.7762638          1.5 drug
## 5  0.236 -0.18892048  0.5897263 0.236 -0.19613539  0.5917346 0.382 -0.48798116   0.7435092              drug
## 6  0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192      -0.5 neutral
## 7  0.312 -0.10857039  0.6407610 0.357 -0.06553960  0.6714658 0.527 -0.14027262   0.8034454      -1.5 neutral
## 8  0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192       0.5 neutral
## 9  0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192       1.5 neutral
## 10 0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192           neutral
## 11 0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192  -0.5DrugvNeutral
## 12 0.508  0.12668409  0.7598382 0.508  0.11930479  0.7611390 0.673  0.21317660   0.8643713  -1.5DrugvNeutral
## 13 0.177 -0.24786607  0.5478653 0.177 -0.25488525  0.5500205 0.301 -0.68415032   0.7096945   0.5DrugvNeutral
## 14 0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192   1.5DrugvNeutral
## 15 0.044 -0.36929027  0.4469343 0.044 -0.37574316  0.4493995 0.085 -1.20380949   0.6201182 DrugvNeutralEarly
## 16 0.385 -0.02560338  0.6873182 0.385 -0.03308957  0.6889423 0.556 -0.06844391   0.8158269  DrugvNeutralLate
## 17 0.305 -0.11660435  0.6359398 0.305 -0.12398856  0.6377734 0.467 -0.28307521   0.7788298      DrugvNeutral
## Saving 7 x 5 in image

## [1] "###tDCS MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 267.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.4414 -0.5660  0.0246  0.5511  4.0628 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.002227 0.04719 
##  id           (Intercept) 0.018720 0.13682 
##  Residual                 0.081454 0.28540 
## Number of obs: 520, groups:  condition:id, 130; id, 65
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)         -0.004655   0.047007  73.685426  -0.099   0.9214    
## conditiondrug       -0.147297   0.026365  63.999993  -5.587 5.08e-07 ***
## time                -0.037129   0.015831 387.999987  -2.345   0.0195 *  
## motion              -0.570243   0.387545  62.999998  -1.471   0.1462    
## conditiondrug:time  -0.009138   0.022389 387.999987  -0.408   0.6834    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.280                     
## time         0.000  0.000              
## motion      -0.844  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000
## [1] "###Neurocaps MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 264.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.3714 -0.3807  0.0382  0.3826  4.7089 
## 
## Random effects:
##  Groups       Name        Variance  Std.Dev. 
##  condition:id (Intercept) 8.537e-17 9.239e-09
##  id           (Intercept) 3.319e-02 1.822e-01
##  Residual                 1.536e-01 3.919e-01
## Number of obs: 232, groups:  condition:id, 58; id, 29
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)         -0.07436    0.09878  31.05170  -0.753  0.45725   
## conditiondrug       -0.16547    0.05146 200.00000  -3.215  0.00152 **
## time                -0.01313    0.03255 200.00000  -0.403  0.68713   
## motion              -0.03058    0.86533  27.00000  -0.035  0.97207   
## conditiondrug:time  -0.01174    0.04603 200.00000  -0.255  0.79902   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.260                     
## time         0.000  0.000              
## motion      -0.864  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000
## [1] "###Neurocaps OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 140.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8340 -0.3946  0.0703  0.3986  4.0753 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.02182  0.1477  
##  id           (Intercept) 0.01711  0.1308  
##  Residual                 0.09622  0.3102  
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)  
## (Intercept)         -0.05875    0.09447  25.22136  -0.622   0.5396  
## conditiondrug       -0.14793    0.06458  21.00000  -2.291   0.0324 *
## time                -0.02566    0.02958 130.00000  -0.868   0.3871  
## motion               0.29626    0.72991  20.00000   0.406   0.6891  
## conditiondrug:time  -0.01035    0.04183 130.00000  -0.247   0.8050  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.342                     
## time         0.000  0.000              
## motion      -0.824  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000
## [1] "###tACS OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 167.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8894 -0.2979  0.0289  0.4702  3.4795 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.008613 0.09281 
##  id           (Intercept) 0.025896 0.16092 
##  Residual                 0.119896 0.34626 
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)  
## (Intercept)         -0.087081   0.108307  23.218330  -0.804   0.4295  
## conditiondrug       -0.161470   0.059228  20.999998  -2.726   0.0127 *
## time                -0.016560   0.033015 130.000003  -0.502   0.6168  
## motion               0.047493   1.051324  20.000002   0.045   0.9644  
## conditiondrug:time  -0.005481   0.046690 130.000003  -0.117   0.9067  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.273                     
## time         0.000  0.000              
## motion      -0.866  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000

ConditionXTime

VMPFC

knitr::include_graphics("VMPFC_zm11.png")

plot_combination('1', 'VMPFC', limits_y = c(-0.4, 0.5))
## [1] "###Condition z-scores###"
##                     tdcs neurocaps_mcr neurocaps_ocr     tacs
## tdcs                  NA    0.07598867     0.5455601 2.524548
## neurocaps_mcr 0.07598867            NA     0.4319498 2.190228
## neurocaps_ocr 0.54556014    0.43194985            NA 1.565483
## tacs          2.52454835    2.19022832     1.5654825       NA
## [1] "###Condition by time z-scores###"
##                     tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                  NA    -0.7157979     0.4934780 -2.132664
## neurocaps_mcr -0.7157979            NA     0.9771868 -1.340874
## neurocaps_ocr  0.4934780     0.9771868            NA -2.096726
## tacs          -2.1326637    -1.3408743    -2.0967259        NA
## [1] "###Condition p-values###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA    0.93942811     0.5853683 0.01158471
## neurocaps_mcr 0.93942811            NA     0.6657779 0.02850768
## neurocaps_ocr 0.58536833    0.66577786            NA 0.11746981
## tacs          0.01158471    0.02850768     0.1174698         NA
## [1] "###Condition by time p-values###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA     0.4741162    0.62167491 0.03295233
## neurocaps_mcr 0.47411617            NA    0.32847669 0.17996127
## neurocaps_ocr 0.62167491     0.3284767            NA 0.03601785
## tacs          0.03295233     0.1799613    0.03601785         NA
## Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

##      ICC1  ICC1_lower ICC1_upper  ICC3  ICC3_lower ICC3_upper ICC3k ICC3k_lower ICC3k_upper               var
## 1   0.073 -0.34407149  0.4697585 0.073 -0.35066003  0.4721591 0.136 -1.08005068   0.6414511         -0.5 drug
## 2   0.079 -0.33889315  0.4743142 0.223 -0.20974246  0.5824346 0.364 -0.53082051   0.7361247         -1.5 drug
## 3   0.042 -0.37171705  0.4446805 0.042 -0.37815638  0.4471519 0.080 -1.21624271   0.6179750          0.5 drug
## 4   0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192          1.5 drug
## 5   0.114 -0.30707170  0.5013852 0.114 -0.31384216  0.5036910 0.205 -0.91478125   0.6699395              drug
## 6  -0.021 -0.42418129  0.3932883 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192      -0.5 neutral
## 7   0.334 -0.08416782  0.6550517 0.345 -0.08013093  0.6633308 0.512 -0.17422246   0.7975933      -1.5 neutral
## 8  -0.056 -0.45308572  0.3626098 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192       0.5 neutral
## 9   0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192       1.5 neutral
## 10 -0.037 -0.43721342  0.3796780 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192           neutral
## 11 -0.001 -0.40777027  0.4099324 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192  -0.5DrugvNeutral
## 12  0.121 -0.30082855  0.5065169 0.314 -0.11408880  0.6436927 0.478 -0.25756262   0.7832276  -1.5DrugvNeutral
## 13 -0.002 -0.40857976  0.4091240 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192   0.5DrugvNeutral
## 14  0.431  0.02938104  0.7152510 0.448  0.04316295  0.7269512 0.619  0.08275399   0.8418897   1.5DrugvNeutral
## 15 -0.048 -0.44653430  0.3697218 0.070 -0.35369437  0.4694629 0.130 -1.09451116   0.6389585 DrugvNeutralEarly
## 16  0.227 -0.19843003  0.5832447 0.227 -0.20561686  0.5852764 0.370 -0.51767681   0.7383904  DrugvNeutralLate
## 17  0.234 -0.19073916  0.5884945 0.311 -0.11713950  0.6418783 0.475 -0.26536355   0.7818829      DrugvNeutral
## Saving 7 x 5 in image

## [1] "###tDCS MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 760.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.4104 -0.5039 -0.0168  0.5026  5.2248 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.003389 0.05821 
##  id           (Intercept) 0.038847 0.19710 
##  Residual                 0.217516 0.46639 
## Number of obs: 520, groups:  condition:id, 130; id, 65
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          0.05693    0.07098  75.08547   0.802   0.4251    
## conditiondrug        0.19212    0.04216  64.00001   4.557 2.39e-05 ***
## time                 0.05140    0.02587 388.00001   1.987   0.0477 *  
## motion              -0.99171    0.58216  63.00000  -1.703   0.0934 .  
## conditiondrug:time  -0.18186    0.03659 388.00001  -4.971 1.00e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.297                     
## time         0.000  0.000              
## motion      -0.840  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000
## [1] "###Neurocaps MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 315.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0960 -0.5486  0.0200  0.5157  3.4532 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.000841 0.0290  
##  id           (Intercept) 0.013410 0.1158  
##  Residual                 0.205415 0.4532  
## Number of obs: 232, groups:  condition:id, 58; id, 29
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)         -0.04938    0.08809  33.98498  -0.561  0.57873   
## conditiondrug        0.18654    0.06000  27.99995   3.109  0.00428 **
## time                 0.08372    0.03764 171.99995   2.224  0.02743 * 
## motion              -0.02502    0.75146  27.00000  -0.033  0.97369   
## conditiondrug:time  -0.13563    0.05323 171.99995  -2.548  0.01171 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.341                     
## time         0.000  0.000              
## motion      -0.842  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000
## [1] "###Neurocaps OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 278.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0190 -0.5387 -0.0104  0.6393  3.7770 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.00000  0.0000  
##  id           (Intercept) 0.01993  0.1412  
##  Residual                 0.25327  0.5033  
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)          0.06942    0.10767  26.00213   0.645  0.52474   
## conditiondrug        0.14476    0.07587 151.00000   1.908  0.05828 . 
## time                 0.06562    0.04798 151.00000   1.367  0.17351   
## motion              -1.40506    0.82846  20.00000  -1.696  0.10540   
## conditiondrug:time  -0.21991    0.06786 151.00000  -3.241  0.00147 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.352                     
## time         0.000  0.000              
## motion      -0.821  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000
## [1] "###tACS OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 266.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1246 -0.5214 -0.0039  0.5059  4.4048 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.0000   0.0000  
##  id           (Intercept) 0.0325   0.1803  
##  Residual                 0.2307   0.4804  
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)          0.078451   0.126663  23.693024   0.619    0.542
## conditiondrug       -0.019427   0.072416 151.000000  -0.268    0.789
## time                 0.006036   0.045800 151.000000   0.132    0.895
## motion              -0.274236   1.224876  20.000000  -0.224    0.825
## conditiondrug:time  -0.023216   0.064770 151.000000  -0.358    0.721
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.286                     
## time         0.000  0.000              
## motion      -0.863  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000

LSTG

knitr::include_graphics("LSFG_z11.png")

plot_combination('3', 'LSTG', limits_y = c(-0.5, 0.3))
## [1] "###Condition z-scores###"
##                   tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                NA     1.9347464     1.0097603  1.2324662
## neurocaps_mcr 1.934746            NA    -0.7885670 -0.6726636
## neurocaps_ocr 1.009760    -0.7885670            NA  0.1405627
## tacs          1.232466    -0.6726636     0.1405627         NA
## [1] "###Condition by time z-scores###"
##                     tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                  NA    -0.7002844     0.3893853 -2.225859
## neurocaps_mcr -0.7002844            NA     0.8720940 -1.115811
## neurocaps_ocr  0.3893853     0.8720940            NA -2.073894
## tacs          -2.2258591    -1.1158114    -2.0738943        NA
## [1] "###Condition p-values###"
##                     tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                  NA    0.05302142     0.3126101 0.2177750
## neurocaps_mcr 0.05302142            NA     0.4303651 0.5011613
## neurocaps_ocr 0.31261014    0.43036512            NA 0.8882154
## tacs          0.21777498    0.50116133     0.8882154        NA
## [1] "###Condition by time p-values###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA     0.4837497    0.69699112 0.02602362
## neurocaps_mcr 0.48374969            NA    0.38315709 0.26450288
## neurocaps_ocr 0.69699112     0.3831571            NA 0.03808913
## tacs          0.02602362     0.2645029    0.03808913         NA
## Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

##      ICC1 ICC1_lower ICC1_upper  ICC3 ICC3_lower ICC3_upper ICC3k ICC3k_lower ICC3k_upper               var
## 1   0.000 -0.4070147  0.4106858 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192         -0.5 drug
## 2   0.122 -0.3002244  0.5070105 0.146 -0.2846024  0.5272714 0.254  -0.7956483   0.6904750         -1.5 drug
## 3   0.215 -0.2106148  0.5747915 0.219 -0.2138283  0.5796010 0.359  -0.5439736   0.7338575          0.5 drug
## 4   0.000 -0.4070147  0.4106858 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192          1.5 drug
## 5   0.078 -0.3394387  0.4738363 0.078 -0.3460512  0.4762250 0.145  -1.0583433   0.6451929              drug
## 6   0.000 -0.4070147  0.4106858 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192      -0.5 neutral
## 7   0.135 -0.2874377  0.5173342 0.217 -0.2156613  0.5783235 0.357  -0.5499188   0.7328326      -1.5 neutral
## 8  -0.002 -0.4084223  0.4092813 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192       0.5 neutral
## 9   0.038 -0.3747622  0.4418376 0.055 -0.3664769  0.4579340 0.104  -1.1569490   0.6281958       1.5 neutral
## 10  0.000 -0.4070147  0.4106858 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192           neutral
## 11  0.000 -0.4070147  0.4106858 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192  -0.5DrugvNeutral
## 12 -0.009 -0.4148260  0.4028425 0.126 -0.3030504  0.5125382 0.224  -0.8696480   0.6777193  -1.5DrugvNeutral
## 13  0.113 -0.3082963  0.5003719 0.113 -0.3150611  0.5026808 0.203  -0.9199683   0.6690453   0.5DrugvNeutral
## 14 -0.003 -0.4096944  0.4080088 0.007 -0.4070884  0.4193681 0.015  -1.3731840   0.5909223   1.5DrugvNeutral
## 15 -0.009 -0.4142570  0.4034180 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192 DrugvNeutralEarly
## 16  0.110 -0.3107377  0.4983452 0.110 -0.3174910  0.5006603 0.199  -0.9303644   0.6672533  DrugvNeutralLate
## 17  0.282 -0.1409204  0.6209804 0.282 -0.1482564  0.6228720 0.440  -0.3481246   0.7676169      DrugvNeutral
## Saving 7 x 5 in image

## [1] "###tDCS MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 488
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.7012 -0.5367  0.0335  0.5085  4.7125 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.0000   0.0000  
##  id           (Intercept) 0.0166   0.1288  
##  Residual                 0.1326   0.3642  
## Number of obs: 520, groups:  condition:id, 130; id, 65
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)         -0.008367   0.050056  77.948844  -0.167   0.8677    
## conditiondrug       -0.043988   0.031941 452.000000  -1.377   0.1691    
## time                 0.100070   0.020201 452.000000   4.954 1.03e-06 ***
## motion              -0.801334   0.407467  63.000000  -1.967   0.0536 .  
## conditiondrug:time  -0.152424   0.028569 452.000000  -5.335 1.51e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.319                     
## time         0.000  0.000              
## motion      -0.833  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000
## [1] "###Neurocaps MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 352.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.4520 -0.4321 -0.0093  0.3784  8.6760 
## 
## Random effects:
##  Groups       Name        Variance  Std.Dev. 
##  condition:id (Intercept) 9.759e-16 3.124e-08
##  id           (Intercept) 2.711e-02 1.646e-01
##  Residual                 2.358e-01 4.856e-01
## Number of obs: 232, groups:  condition:id, 58; id, 29
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)         -0.16974    0.10412  32.82824  -1.630  0.11262   
## conditiondrug       -0.18197    0.06376 200.00000  -2.854  0.00477 **
## time                 0.07866    0.04033 200.00000   1.951  0.05250 . 
## motion               1.86176    0.89932  27.00000   2.070  0.04812 * 
## conditiondrug:time  -0.10776    0.05703 200.00000  -1.889  0.06028 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.306                     
## time         0.000  0.000              
## motion      -0.852  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000
## [1] "###Neurocaps OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 204.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7379 -0.4146  0.0566  0.5143  3.8043 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.00000  0.0000  
##  id           (Intercept) 0.02198  0.1482  
##  Residual                 0.15951  0.3994  
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)          0.007111   0.095683  24.596763   0.074 0.941362    
## conditiondrug       -0.112810   0.060209 151.000000  -1.874 0.062913 .  
## time                 0.131431   0.038080 151.000000   3.451 0.000723 ***
## motion              -0.286169   0.746753  20.000000  -0.383 0.705601    
## conditiondrug:time  -0.176162   0.053853 151.000000  -3.271 0.001327 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.315                     
## time         0.000  0.000              
## motion      -0.832  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000
## [1] "###tACS OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 201.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.2358 -0.4387  0.0239  0.4450  3.5822 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.00000  0.0000  
##  id           (Intercept) 0.07051  0.2655  
##  Residual                 0.14269  0.3777  
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)  
## (Intercept)         -0.03620    0.14842  21.55345  -0.244   0.8096  
## conditiondrug       -0.12446    0.05695 151.00000  -2.186   0.0304 *
## time                 0.01476    0.03602 151.00000   0.410   0.6825  
## motion               0.31504    1.46995  20.00000   0.214   0.8325  
## conditiondrug:time  -0.02244    0.05093 151.00000  -0.440   0.6602  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.192                     
## time         0.000  0.000              
## motion      -0.884  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000

RSTG

knitr::include_graphics("RSFG_z7.png")

plot_combination('4', 'RSTG', limits_y = c(-0.4, 0.25))
## [1] "###Condition z-scores###"
##                    tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                 NA     2.2100759     0.3935193  1.6060658
## neurocaps_mcr 2.2100759            NA    -1.5704733 -0.2390023
## neurocaps_ocr 0.3935193    -1.5704733            NA  1.1447217
## tacs          1.6060658    -0.2390023     1.1447217         NA
## [1] "###Condition by time z-scores###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA    -1.3968480    -0.6932403 -2.4421241
## neurocaps_mcr -1.3968480            NA     0.5854781 -0.9980761
## neurocaps_ocr -0.6932403     0.5854781            NA -1.5448472
## tacs          -2.4421241    -0.9980761    -1.5448472         NA
## [1] "###Condition p-values###"
##                    tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                 NA     0.0270999     0.6939360 0.1082595
## neurocaps_mcr 0.0270999            NA     0.1163050 0.8111038
## neurocaps_ocr 0.6939360     0.1163050            NA 0.2523245
## tacs          0.1082595     0.8111038     0.2523245        NA
## [1] "###Condition by time p-values###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA     0.1624593     0.4881588 0.01460113
## neurocaps_mcr 0.16245930            NA     0.5582263 0.31824245
## neurocaps_ocr 0.48815879     0.5582263            NA 0.12238323
## tacs          0.01460113     0.3182424     0.1223832         NA
## Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

##      ICC1  ICC1_lower ICC1_upper  ICC3  ICC3_lower ICC3_upper ICC3k ICC3k_lower ICC3k_upper               var
## 1   0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192         -0.5 drug
## 2   0.395 -0.01412295  0.6933300 0.395 -0.02161320  0.6949285 0.566 -0.04418130   0.8200092         -1.5 drug
## 3   0.386 -0.02518190  0.6875407 0.386 -0.03266827  0.6891638 0.557 -0.06754306   0.8159822          0.5 drug
## 4   0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192          1.5 drug
## 5   0.165 -0.25966198  0.5389808 0.165 -0.26663577  0.5411658 0.283 -0.72715781   0.7022811              drug
## 6   0.067 -0.34948960  0.4649449 0.067 -0.35604980  0.4673594 0.126 -1.10583024   0.6370074      -0.5 neutral
## 7   0.070 -0.34684083  0.4673042 0.090 -0.33541924  0.4854722 0.166 -1.00941603   0.6536268      -1.5 neutral
## 8   0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192       0.5 neutral
## 9   0.034 -0.37789252  0.4388980 0.036 -0.38257773  0.4430066 0.070 -1.23927418   0.6140050       1.5 neutral
## 10  0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192           neutral
## 11 -0.012 -0.41688235  0.4007577 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192  -0.5DrugvNeutral
## 12  0.116 -0.30546157  0.5027142 0.166 -0.26578608  0.5418121 0.284 -0.72400175   0.7028251  -1.5DrugvNeutral
## 13  0.193 -0.23183343  0.5596598 0.193 -0.23891098  0.5617747 0.324 -0.62781349   0.7194056   0.5DrugvNeutral
## 14  0.050 -0.36450590  0.4513477 0.050 -0.37098524  0.4538006 0.095 -1.17957565   0.6242955   1.5DrugvNeutral
## 15 -0.003 -0.40969686  0.4080063 0.050 -0.37129867  0.4535119 0.094 -1.18116075   0.6240223 DrugvNeutralEarly
## 16  0.167 -0.25744436  0.5406647 0.167 -0.26442684  0.5428441 0.286 -0.71896816   0.7036928  DrugvNeutralLate
## 17  0.283 -0.13993787  0.6215959 0.304 -0.12516779  0.6370622 0.466 -0.28615267   0.7782994      DrugvNeutral
## Saving 7 x 5 in image

## [1] "###tDCS MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 460.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.2027 -0.4741  0.0263  0.5178  4.4562 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.00000  0.0000  
##  id           (Intercept) 0.01976  0.1406  
##  Residual                 0.12372  0.3517  
## Number of obs: 520, groups:  condition:id, 130; id, 65
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)          0.006044   0.051259  76.064818   0.118    0.906    
## conditiondrug       -0.037236   0.030849 451.999997  -1.207    0.228    
## time                 0.094633   0.019511 451.999997   4.850 1.70e-06 ***
## motion              -0.530157   0.419862  63.000004  -1.263    0.211    
## conditiondrug:time  -0.152584   0.027593 451.999997  -5.530 5.43e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.301                     
## time         0.000  0.000              
## motion      -0.839  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000
## [1] "###Neurocaps MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 238.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5813 -0.4029  0.0069  0.4913  5.5167 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.006746 0.08213 
##  id           (Intercept) 0.028415 0.16857 
##  Residual                 0.132720 0.36431 
## Number of obs: 232, groups:  condition:id, 58; id, 29
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)         -0.13853    0.09533  31.40335  -1.453  0.15610   
## conditiondrug       -0.17176    0.05247  28.00002  -3.273  0.00283 **
## time                 0.05542    0.03025 172.00005   1.832  0.06869 . 
## motion               1.17457    0.83155  27.00000   1.413  0.16923   
## conditiondrug:time  -0.08147    0.04279 172.00005  -1.904  0.05857 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.275                     
## time         0.000  0.000              
## motion      -0.861  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000
## [1] "###Neurocaps OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 129.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8687 -0.4822  0.0144  0.6112  4.2070 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.00000  0.0000  
##  id           (Intercept) 0.01532  0.1238  
##  Residual                 0.10256  0.3202  
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)          0.06867    0.07824  24.39391   0.878  0.38869   
## conditiondrug       -0.05978    0.04828 151.00000  -1.238  0.21754   
## time                 0.09106    0.03053 151.00000   2.982  0.00334 **
## motion              -0.69649    0.61190  20.00000  -1.138  0.26847   
## conditiondrug:time  -0.11706    0.04318 151.00000  -2.711  0.00749 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.309                     
## time         0.000  0.000              
## motion      -0.834  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000
## [1] "###tACS OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 187.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1370 -0.4503  0.0529  0.4722  4.2635 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.01385  0.1177  
##  id           (Intercept) 0.04863  0.2205  
##  Residual                 0.12704  0.3564  
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)  
## (Intercept)         -0.030549   0.134884  22.411828  -0.226   0.8229  
## conditiondrug       -0.151911   0.064393  20.999999  -2.359   0.0281 *
## time                -0.007535   0.033984 129.999999  -0.222   0.8249  
## motion               0.688983   1.321838  20.000001   0.521   0.6079  
## conditiondrug:time  -0.017246   0.048061 129.999999  -0.359   0.7203  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.239                     
## time         0.000  0.000              
## motion      -0.874  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000

LVStriatum

knitr::include_graphics("LVStri_zm1.png")

plot_combination('5', 'LVSTRI', limits_y = c(-0.22, 0.35))
## [1] "###Condition z-scores###"
##                    tdcs neurocaps_mcr neurocaps_ocr     tacs
## tdcs                 NA     0.7779924     0.1985035 2.713589
## neurocaps_mcr 0.7779924            NA    -0.3694164 1.641814
## neurocaps_ocr 0.1985035    -0.3694164            NA 1.762054
## tacs          2.7135890     1.6418139     1.7620536       NA
## [1] "###Condition by time z-scores###"
##                     tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                  NA     0.1731414     0.9808995 -2.384295
## neurocaps_mcr  0.1731414            NA     0.7411386 -2.132148
## neurocaps_ocr  0.9808995     0.7411386            NA -2.544832
## tacs          -2.3842952    -2.1321485    -2.5448317        NA
## [1] "###Condition p-values###"
##                      tdcs neurocaps_mcr neurocaps_ocr        tacs
## tdcs                   NA     0.4365735    0.84265117 0.006655868
## neurocaps_mcr 0.436573519            NA    0.71181736 0.100628576
## neurocaps_ocr 0.842651170     0.7118174            NA 0.078060239
## tacs          0.006655868     0.1006286    0.07806024          NA
## [1] "###Condition by time p-values###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA    0.86254031    0.32664230 0.01711187
## neurocaps_mcr 0.86254031            NA    0.45860940 0.03299464
## neurocaps_ocr 0.32664230    0.45860940            NA 0.01093304
## tacs          0.01711187    0.03299464    0.01093304         NA
## Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

##      ICC1 ICC1_lower ICC1_upper  ICC3  ICC3_lower ICC3_upper ICC3k ICC3k_lower ICC3k_upper               var
## 1  -0.011 -0.4163157  0.4013330 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192         -0.5 drug
## 2  -0.233 -0.5842313  0.1969920 0.048 -0.37265010  0.4522652 0.092 -1.18801360   0.6228410         -1.5 drug
## 3   0.180 -0.2452362  0.5498219 0.180 -0.25226520  0.5519705 0.305 -0.67474510   0.7113157          0.5 drug
## 4   0.000 -0.4070147  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192          1.5 drug
## 5  -0.042 -0.4417737  0.3748305 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192              drug
## 6   0.281 -0.1420422  0.6202767 0.281 -0.14937578  0.6221709 0.439 -0.35121451   0.7670843      -0.5 neutral
## 7   0.128 -0.2940845  0.5119968 0.129 -0.30031985  0.5147496 0.228 -0.85844898   0.6796497      -1.5 neutral
## 8   0.000 -0.4070147  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192       0.5 neutral
## 9   0.000 -0.4070147  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192       1.5 neutral
## 10  0.000 -0.4070147  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192           neutral
## 11  0.000 -0.4070147  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192  -0.5DrugvNeutral
## 12  0.050 -0.3645727  0.4512863 0.440  0.03306215  0.7221461 0.611  0.06400805   0.8386583  -1.5DrugvNeutral
## 13  0.165 -0.2592465  0.5392968 0.165 -0.26622191  0.5414807 0.284 -0.72561967   0.7025462   0.5DrugvNeutral
## 14  0.202 -0.2236710  0.5655438 0.202 -0.23077683  0.5676383 0.336 -0.60002570   0.7241955   1.5DrugvNeutral
## 15 -0.131 -0.5104297  0.2960213 0.125 -0.30351702  0.5121592 0.223 -0.87157053   0.6773879 DrugvNeutralEarly
## 16  0.236 -0.1889639  0.5896970 0.236 -0.19617868  0.5917054 0.382 -0.48811515   0.7434861  DrugvNeutralLate
## 17  0.218 -0.2078345  0.5767352 0.317 -0.11044641  0.6458480 0.482 -0.24831874   0.7848210      DrugvNeutral
## Saving 7 x 5 in image

## [1] "###tDCS MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 153.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5172 -0.5441 -0.0123  0.5962  4.8206 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.000000 0.00000 
##  id           (Intercept) 0.008348 0.09137 
##  Residual                 0.069458 0.26355 
## Number of obs: 520, groups:  condition:id, 130; id, 65
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          0.09502    0.03590  78.26601   2.647 0.009815 ** 
## conditiondrug        0.11960    0.02311 452.00000   5.174 3.45e-07 ***
## time                 0.02737    0.01462 452.00000   1.872 0.061827 .  
## motion              -1.21634    0.29192  63.00000  -4.167 9.59e-05 ***
## conditiondrug:time  -0.07181    0.02067 452.00000  -3.473 0.000564 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.322                     
## time         0.000  0.000              
## motion      -0.832  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000
## [1] "###Neurocaps MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 137.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.7549 -0.4475  0.0024  0.4298  6.3774 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.00000  0.0000  
##  id           (Intercept) 0.02711  0.1647  
##  Residual                 0.08492  0.2914  
## Number of obs: 232, groups:  condition:id, 58; id, 29
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)  
## (Intercept)         -0.13960    0.08317  30.08911  -1.679   0.1036  
## conditiondrug        0.08482    0.03826 200.00000   2.217   0.0278 *
## time                 0.05273    0.02420 200.00000   2.179   0.0305 *
## motion               0.95340    0.73434  27.00000   1.298   0.2052  
## conditiondrug:time  -0.07873    0.03422 200.00000  -2.300   0.0225 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.230                     
## time         0.000  0.000              
## motion      -0.871  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000
## [1] "###Neurocaps OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 154.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.4172 -0.5649  0.0638  0.5921  3.0391 
## 
## Random effects:
##  Groups       Name        Variance  Std.Dev. 
##  condition:id (Intercept) 9.202e-19 9.592e-10
##  id           (Intercept) 2.537e-02 1.593e-01
##  Residual                 1.155e-01 3.398e-01
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)          0.05760    0.09215  23.46867   0.625  0.53792   
## conditiondrug        0.10844    0.05123 151.00000   2.117  0.03592 * 
## time                 0.03162    0.03240 151.00000   0.976  0.33072   
## motion              -0.72282    0.72777  20.00000  -0.993  0.33248   
## conditiondrug:time  -0.12112    0.04582 151.00000  -2.643  0.00908 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.278                     
## time         0.000  0.000              
## motion      -0.842  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000
## [1] "###tACS OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 73.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.3930 -0.5325  0.0113  0.5596  2.3392 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.00000  0.0000  
##  id           (Intercept) 0.01882  0.1372  
##  Residual                 0.07131  0.2670  
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)         -0.063214   0.084064  22.495706  -0.752    0.460
## conditiondrug       -0.006365   0.040257 150.999999  -0.158    0.875
## time                -0.011709   0.025461 150.999999  -0.460    0.646
## motion               0.540027   0.823652  20.000001   0.656    0.520
## conditiondrug:time   0.027189   0.036007 150.999999   0.755    0.451
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.239                     
## time         0.000  0.000              
## motion      -0.874  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000

RVStriatum

knitr::include_graphics("RVStri_zm1.png")

plot_combination('6', 'RVSTRI', limits_y = c(-0.2, 0.3))
## [1] "###Condition z-scores###"
##                    tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                 NA     0.5332700     0.7414731 1.9582658
## neurocaps_mcr 0.5332700            NA     0.2394436 1.2620420
## neurocaps_ocr 0.7414731     0.2394436            NA 0.9295674
## tacs          1.9582658     1.2620420     0.9295674        NA
## [1] "###Condition by time z-scores###"
##                     tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                  NA    -0.7847964     0.7936411 -2.304104
## neurocaps_mcr -0.7847964            NA     1.2817258 -1.361958
## neurocaps_ocr  0.7936411     1.2817258            NA -2.460801
## tacs          -2.3041039    -1.3619583    -2.4608006        NA
## [1] "###Condition p-values###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA     0.5938467     0.4584066 0.05019883
## neurocaps_mcr 0.59384666            NA     0.8107616 0.20693369
## neurocaps_ocr 0.45840665     0.8107616            NA 0.35259512
## tacs          0.05019883     0.2069337     0.3525951         NA
## [1] "###Condition by time p-values###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA     0.4325729    0.42740442 0.02121681
## neurocaps_mcr 0.43257295            NA    0.19993887 0.17321106
## neurocaps_ocr 0.42740442     0.1999389            NA 0.01386274
## tacs          0.02121681     0.1732111    0.01386274         NA
## Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

##      ICC1  ICC1_lower ICC1_upper  ICC3  ICC3_lower ICC3_upper ICC3k ICC3k_lower ICC3k_upper               var
## 1  -0.013 -0.41814081  0.3994775 0.000 -0.41324703  0.4132470 0.000  -1.4085895   0.5848192         -0.5 drug
## 2  -0.183 -0.54883570  0.2465632 0.000 -0.41324703  0.4132470 0.000  -1.4085895   0.5848192         -1.5 drug
## 3   0.173 -0.25189290  0.5448524 0.173 -0.25889682  0.5470178 0.295  -0.6986795   0.7071900          0.5 drug
## 4   0.209 -0.21594380  0.5710411 0.209 -0.22307541  0.5731163 0.346  -0.5742524   0.7286381          1.5 drug
## 5  -0.004 -0.41043615  0.4072652 0.000 -0.41324703  0.4132470 0.000  -1.4085895   0.5848192              drug
## 6   0.466  0.07263964  0.7357876 0.466  0.06518258  0.7371991 0.636   0.1223876   0.8487215      -0.5 neutral
## 7   0.183 -0.24181183  0.5523566 0.183 -0.24885350  0.5544966 0.310  -0.6625964   0.7134098      -1.5 neutral
## 8   0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000  -1.4085895   0.5848192       0.5 neutral
## 9   0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000  -1.4085895   0.5848192       1.5 neutral
## 10  0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000  -1.4085895   0.5848192           neutral
## 11 -0.007 -0.41276396  0.4049248 0.000 -0.41324703  0.4132470 0.000  -1.4085895   0.5848192  -0.5DrugvNeutral
## 12 -0.103 -0.48941402  0.3213702 0.030 -0.38776739  0.4380956 0.059  -1.2667322   0.6092719  -1.5DrugvNeutral
## 13  0.204 -0.22138222  0.5671794 0.204 -0.22849574  0.5692682 0.339  -0.5923383   0.7255206   0.5DrugvNeutral
## 14  0.081 -0.33694092  0.4760205 0.093 -0.33326886  0.4873205 0.170  -0.9997099   0.6552999   1.5DrugvNeutral
## 15 -0.089 -0.47856170  0.3340204 0.046 -0.37437359  0.4506705 0.088  -1.1967960   0.6213272 DrugvNeutralEarly
## 16  0.149 -0.27471736  0.5273771 0.149 -0.28163020  0.5296003 0.259  -0.7840814   0.6924689  DrugvNeutralLate
## 17  0.308 -0.11316548  0.6380107 0.319 -0.10837702  0.6470672 0.484  -0.2431006   0.7857205      DrugvNeutral
## Saving 7 x 5 in image

## [1] "###tDCS MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 80.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.9640 -0.4675  0.0041  0.5317  5.0395 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.00000  0.00000 
##  id           (Intercept) 0.00823  0.09072 
##  Residual                 0.05978  0.24450 
## Number of obs: 520, groups:  condition:id, 130; id, 65
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          0.05576    0.03435  77.20889   1.623 0.108602    
## conditiondrug        0.10634    0.02144 452.00000   4.959    1e-06 ***
## time                 0.01734    0.01356 452.00000   1.279 0.201721    
## motion              -0.99538    0.28031  63.00000  -3.551 0.000733 ***
## conditiondrug:time  -0.06362    0.01918 452.00000  -3.317 0.000984 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.312                     
## time         0.000  0.000              
## motion      -0.835  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000
## [1] "###Neurocaps MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 75.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7190 -0.5456 -0.0613  0.5309  5.0944 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.00000  0.0000  
##  id           (Intercept) 0.01216  0.1103  
##  Residual                 0.06758  0.2600  
## Number of obs: 232, groups:  condition:id, 58; id, 29
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)  
## (Intercept)         -0.14029    0.06220  31.54753  -2.255   0.0312 *
## conditiondrug        0.08485    0.03413 200.00000   2.486   0.0137 *
## time                 0.01649    0.02159 200.00000   0.764   0.4458  
## motion               0.61150    0.54268  27.00000   1.127   0.2697  
## conditiondrug:time  -0.03532    0.03053 200.00000  -1.157   0.2487  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.274                     
## time         0.000  0.000              
## motion      -0.861  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000
## [1] "###Neurocaps OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 76.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3723 -0.5367  0.0829  0.5963  2.4074 
## 
## Random effects:
##  Groups       Name        Variance  Std.Dev. 
##  condition:id (Intercept) 1.243e-15 3.526e-08
##  id           (Intercept) 1.485e-02 1.218e-01
##  Residual                 7.369e-02 2.715e-01
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)          0.03235    0.07179  23.66917   0.451  0.65636   
## conditiondrug        0.07209    0.04092 151.00000   1.762  0.08017 . 
## time                 0.00827    0.02588 151.00000   0.320  0.74977   
## motion              -0.56711    0.56573  20.00000  -1.002  0.32811   
## conditiondrug:time  -0.09641    0.03660 151.00000  -2.634  0.00932 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.285                     
## time         0.000  0.000              
## motion      -0.841  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000
## [1] "###tACS OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 61
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.3285 -0.4517  0.0488  0.5717  2.5342 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.00000  0.0000  
##  id           (Intercept) 0.01882  0.1372  
##  Residual                 0.06581  0.2565  
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)  
## (Intercept)         -0.15060    0.08288  22.35899  -1.817   0.0826 .
## conditiondrug        0.01975    0.03867 151.00000   0.511   0.6104  
## time                -0.02562    0.02446 151.00000  -1.047   0.2966  
## motion               1.20248    0.81333  20.00000   1.478   0.1549  
## conditiondrug:time   0.02752    0.03459 151.00000   0.796   0.4276  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.233                     
## time         0.000  0.000              
## motion      -0.876  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000

RAmy

knitr::include_graphics("RAmy_zm16.png")

plot_combination('8', 'RAmy', limits_y = c(-0.05, 0.75))
## [1] "###Condition z-scores###"
##                    tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                 NA     0.7234911     2.4657736 2.7482396
## neurocaps_mcr 0.7234911            NA     1.3756300 1.8963280
## neurocaps_ocr 2.4657736     1.3756300            NA 0.7971716
## tacs          2.7482396     1.8963280     0.7971716        NA
## [1] "###Condition by time z-scores###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA    -0.6276469     0.3958587 -1.4661307
## neurocaps_mcr -0.6276469            NA     0.8888891 -0.8061678
## neurocaps_ocr  0.3958587     0.8888891            NA -1.6096743
## tacs          -1.4661307    -0.8061678    -1.6096743         NA
## [1] "###Condition p-values###"
##                      tdcs neurocaps_mcr neurocaps_ocr        tacs
## tdcs                   NA    0.46937821    0.01367177 0.005991622
## neurocaps_mcr 0.469378211            NA    0.16893622 0.057916681
## neurocaps_ocr 0.013671771    0.16893622            NA 0.425351385
## tacs          0.005991622    0.05791668    0.42535139          NA
## [1] "###Condition by time p-values###"
##                    tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                 NA     0.5302353     0.6922093 0.1426127
## neurocaps_mcr 0.5302353            NA     0.3740627 0.4201461
## neurocaps_ocr 0.6922093     0.3740627            NA 0.1074690
## tacs          0.1426127     0.4201461     0.1074690        NA
## Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

##      ICC1 ICC1_lower ICC1_upper  ICC3 ICC3_lower ICC3_upper ICC3k ICC3k_lower ICC3k_upper               var
## 1   0.000 -0.4070147  0.4106858 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192         -0.5 drug
## 2  -0.028 -0.4299547  0.3873026 0.001 -0.4123668  0.4141265 0.002  -1.4034836   0.5856994         -1.5 drug
## 3   0.080 -0.3379906  0.4751038 0.081 -0.3436572  0.4783232 0.150  -1.0471881   0.6471158          0.5 drug
## 4   0.029 -0.3828382  0.4342177 0.029 -0.3892143  0.4367175 0.056  -1.2744711   0.6079379          1.5 drug
## 5   0.000 -0.4070147  0.4106858 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192              drug
## 6  -0.169 -0.5392473  0.2593116 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192      -0.5 neutral
## 7   0.213 -0.2123710  0.5735592 0.213 -0.2195142  0.5756255 0.351  -0.5625066   0.7306628      -1.5 neutral
## 8   0.221 -0.2048928  0.5787821 0.227 -0.2057729  0.5851693 0.370  -0.5181713   0.7383051       0.5 neutral
## 9   0.000 -0.4070147  0.4106858 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192       1.5 neutral
## 10 -0.023 -0.4256577  0.3917642 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192           neutral
## 11  0.000 -0.4070147  0.4106858 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192  -0.5DrugvNeutral
## 12  0.169 -0.2555671  0.5420853 0.214 -0.2181521  0.5765813 0.353  -0.5580423   0.7314324  -1.5DrugvNeutral
## 13  0.295 -0.1275845  0.6292541 0.295 -0.1349480  0.6311138 0.455  -0.3119998   0.7738440   0.5DrugvNeutral
## 14  0.000 -0.4070147  0.4106858 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192   1.5DrugvNeutral
## 15  0.244 -0.1812008  0.5949153 0.324 -0.1030006  0.6502167 0.489  -0.2296560   0.7880380 DrugvNeutralEarly
## 16  0.137 -0.2863589  0.5181947 0.137 -0.2932223  0.5204475 0.240  -0.8297439   0.6845978  DrugvNeutralLate
## 17  0.144 -0.2793620  0.5237359 0.144 -0.2862554  0.5259709 0.252  -0.8021227   0.6893590      DrugvNeutral
## Saving 7 x 5 in image

## [1] "###tDCS MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 596.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3756 -0.5396  0.0180  0.5297  5.3356 
## 
## Random effects:
##  Groups       Name        Variance  Std.Dev.
##  condition:id (Intercept) 0.0007294 0.02701 
##  id           (Intercept) 0.0467436 0.21620 
##  Residual                 0.1527508 0.39083 
## Number of obs: 520, groups:  condition:id, 130; id, 65
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)          0.223221   0.069208  71.367835   3.225 0.001900 ** 
## conditiondrug        0.292493   0.034604  63.999880   8.453 5.15e-12 ***
## time                -0.004978   0.021680 387.999871  -0.230 0.818492    
## motion              -0.568463   0.575556  63.000002  -0.988 0.327090    
## conditiondrug:time  -0.105186   0.030659 387.999871  -3.431 0.000667 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.250                     
## time         0.000  0.000              
## motion      -0.851  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000
## [1] "###Neurocaps MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 219.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6537 -0.6263 -0.0258  0.5302  3.8533 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.01345  0.1160  
##  id           (Intercept) 0.01921  0.1386  
##  Residual                 0.11989  0.3462  
## Number of obs: 232, groups:  condition:id, 58; id, 29
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          0.16587    0.08862  32.64059   1.872 0.070230 .  
## conditiondrug        0.24565    0.05472  28.00000   4.489 0.000112 ***
## time                 0.04765    0.02875 172.00000   1.657 0.099304 .  
## motion               0.02155    0.76478  27.00000   0.028 0.977723    
## conditiondrug:time  -0.07322    0.04066 172.00000  -1.801 0.073518 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.309                     
## time         0.000  0.000              
## motion      -0.851  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000
## [1] "###Neurocaps OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 140.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3272 -0.5281  0.0024  0.5381  3.5504 
## 
## Random effects:
##  Groups       Name        Variance  Std.Dev. 
##  condition:id (Intercept) 3.143e-16 1.773e-08
##  id           (Intercept) 2.639e-02 1.625e-01
##  Residual                 1.056e-01 3.249e-01
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)          0.27838    0.09160  23.17985   3.039  0.00580 **
## conditiondrug        0.14462    0.04898 151.00000   2.953  0.00366 **
## time                 0.01792    0.03098 151.00000   0.578  0.56386   
## motion              -1.16955    0.72570  20.00000  -1.612  0.12271   
## conditiondrug:time  -0.12635    0.04381 151.00000  -2.884  0.00450 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.267                     
## time         0.000  0.000              
## motion      -0.845  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000
## [1] "###tACS OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + motion + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 177.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.89617 -0.49440 -0.00314  0.54474  2.35083 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.02513  0.1585  
##  id           (Intercept) 0.02384  0.1544  
##  Residual                 0.11978  0.3461  
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)  
## (Intercept)          0.24934    0.11658  24.02819   2.139   0.0428 *
## conditiondrug        0.07601    0.07076  21.00000   1.074   0.2949  
## time                -0.04487    0.03300 130.00000  -1.360   0.1762  
## motion              -0.10212    1.12100  20.00000  -0.091   0.9283  
## conditiondrug:time  -0.02332    0.04667 130.00000  -0.500   0.6181  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time   motion
## conditindrg -0.303                     
## time         0.000  0.000              
## motion      -0.858  0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707  0.000

Self-Report and RTs

Craving

plot_combination('craving', 'Craving', limits_y = c(1, 4), imaging = FALSE, forest_range = c(-0.5, 2.5))
## [1] "###Condition z-scores###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA     0.1576549    -3.5832372 -3.6285460
## neurocaps_mcr  0.1576549            NA    -2.8131019 -2.8858501
## neurocaps_ocr -3.5832372    -2.8131019            NA -0.1371176
## tacs          -3.6285460    -2.8858501    -0.1371176         NA
## [1] "###Condition by time z-scores###"
##                     tdcs neurocaps_mcr neurocaps_ocr        tacs
## tdcs                  NA   -1.36784194    -0.6971554 -1.10796210
## neurocaps_mcr -1.3678419            NA     0.4683281  0.08924022
## neurocaps_ocr -0.6971554    0.46832809            NA -0.34567535
## tacs          -1.1079621    0.08924022    -0.3456754          NA
## [1] "###Condition p-values###"
##                       tdcs neurocaps_mcr neurocaps_ocr        tacs
## tdcs                    NA    0.87472870  0.0003393621 0.000285022
## neurocaps_mcr 0.8747287022            NA  0.0049066102 0.003903580
## neurocaps_ocr 0.0003393621    0.00490661            NA 0.890937842
## tacs          0.0002850220    0.00390358  0.8909378416          NA
## [1] "###Condition by time p-values###"
##                    tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                 NA     0.1713616     0.4857055 0.2678782
## neurocaps_mcr 0.1713616            NA     0.6395500 0.9288910
## neurocaps_ocr 0.4857055     0.6395500            NA 0.7295867
## tacs          0.2678782     0.9288910     0.7295867        NA
## Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

##     ICC1  ICC1_lower ICC1_upper  ICC3   ICC3_lower ICC3_upper ICC3k ICC3k_lower ICC3k_upper               var
## 1  0.640  0.31473935  0.8318957 0.648  0.321032941  0.8372305 0.787  0.48603321   0.9114049         -0.5 drug
## 2  0.210 -0.21495602  0.5717387 0.210 -0.222090866  0.5738115 0.348 -0.57099437   0.7291998         -1.5 drug
## 3  0.290 -0.13220264  0.6264082 0.290 -0.139556940  0.6282790 0.450 -0.32438391   0.7717092          0.5 drug
## 4  0.408  0.00172231  0.7014687 0.408 -0.005770302  0.7030322 0.580 -0.01160758   0.8256241          1.5 drug
## 5  0.524  0.14929176  0.7694099 0.524  0.141958008  0.7706655 0.688  0.24862212   0.8704812              drug
## 6  0.593  0.24556160  0.8072723 0.601  0.250377263  0.8126749 0.751  0.40048275   0.8966582      -0.5 neutral
## 7  0.457  0.06109960  0.7304258 0.478  0.081211991  0.7444697 0.647  0.15022399   0.8535198      -1.5 neutral
## 8  0.256 -0.16844405  0.6033506 0.256 -0.175714827  0.6053087 0.408 -0.42634475   0.7541337       0.5 neutral
## 9  0.227 -0.19834601  0.5833024 0.227 -0.205533104  0.5853339 0.370 -0.51741137   0.7384361       1.5 neutral
## 10 0.471  0.07915571  0.7387789 0.471  0.071705701  0.7401768 0.640  0.13381603   0.8506915           neutral
## 11 0.549  0.18239307  0.7829393 0.599  0.246868541  0.8114011 0.749  0.39598167   0.8958823  -0.5DrugvNeutral
## 12 0.169 -0.25556597  0.5420861 0.169 -0.262555756  0.5442607 0.289 -0.71206944   0.7048819  -1.5DrugvNeutral
## 13 0.450  0.05319355  0.7267043 0.456  0.052062876  0.7311302 0.626  0.09897294   0.8446854   0.5DrugvNeutral
## 14 0.585  0.23387444  0.8029099 0.585  0.226779286  0.8040033 0.738  0.36971489   0.8913546   1.5DrugvNeutral
## 15 0.470  0.07798528  0.7382435 0.470  0.070533954  0.7396438 0.640  0.13177341   0.8503394 DrugvNeutralEarly
## 16 0.747  0.48816843  0.8858083 0.753  0.493637658  0.8895781 0.859  0.66098716   0.9415626  DrugvNeutralLate
## 17 0.698  0.40666757  0.8617427 0.698  0.400395025  0.8625346 0.822  0.57183154   0.9261944      DrugvNeutral
## Saving 7 x 5 in image

## [1] "###tDCS MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time +  (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 1014
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.06640 -0.53291  0.05142  0.55519  3.06609 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.1242   0.3524  
##  id           (Intercept) 0.3488   0.5906  
##  Residual                 0.2824   0.5314  
## Number of obs: 494, groups:  condition:id, 129; id, 65
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          1.69379    0.09207  91.37980   18.40  < 2e-16 ***
## conditiondrug        1.31929    0.07862  63.39930   16.78  < 2e-16 ***
## time                 0.18236    0.03106 368.81187    5.87  9.7e-09 ***
## conditiondrug:time  -0.10600    0.04327 367.10038   -2.45   0.0148 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.427              
## time        -0.013  0.013       
## cndtndrg:tm  0.009 -0.015 -0.718
## [1] "###Neurocaps MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time +  (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 370.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3459 -0.4088 -0.0669  0.1513  3.7800 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.36796  0.6066  
##  id           (Intercept) 0.07349  0.2711  
##  Residual                 0.14943  0.3866  
## Number of obs: 228, groups:  condition:id, 58; id, 29
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          1.08915    0.12866  54.96792   8.465 1.54e-11 ***
## conditiondrug        1.29014    0.16737  27.95315   7.708 2.16e-08 ***
## time                 0.04062    0.03221 168.23056   1.261    0.209    
## conditiondrug:time  -0.01991    0.04571 168.28819  -0.436    0.664    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.651              
## time        -0.001  0.001       
## cndtndrg:tm  0.001 -0.003 -0.705
## [1] "###Neurocaps OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time +  (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 296.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8107 -0.5292  0.0609  0.4771  2.8296 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.1905   0.4365  
##  id           (Intercept) 0.1031   0.3212  
##  Residual                 0.1874   0.4329  
## Number of obs: 174, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          1.38243    0.12454  38.34890  11.101 1.53e-13 ***
## conditiondrug        1.91698    0.14711  21.01962  13.031 1.55e-11 ***
## time                 0.14109    0.04182 128.19406   3.374 0.000981 ***
## conditiondrug:time  -0.05492    0.05914 128.18996  -0.929 0.354867    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.591              
## time        -0.007  0.006       
## cndtndrg:tm  0.005 -0.009 -0.707
## [1] "###tACS OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time +  (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 291.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.64887 -0.62209  0.04536  0.61595  2.45138 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.2149   0.4636  
##  id           (Intercept) 0.1002   0.3165  
##  Residual                 0.1777   0.4215  
## Number of obs: 173, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          1.38609    0.12795  39.04699  10.833 2.49e-13 ***
## conditiondrug        1.94616    0.15383  20.95862  12.651 2.81e-11 ***
## time                 0.17284    0.04025 127.16548   4.295 3.44e-05 ***
## conditiondrug:time  -0.02645    0.05730 127.31114  -0.462    0.645    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.600              
## time        -0.002  0.002       
## cndtndrg:tm  0.002 -0.004 -0.702

Craving RT

plot_combination('craving_rt', 'CravingRT', limits_y = c(1, 3.5), imaging = FALSE)
## [1] "###Condition z-scores###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA     -2.385365      1.594898 -0.8158398
## neurocaps_mcr -2.3853647            NA      3.496007  1.4708419
## neurocaps_ocr  1.5948985      3.496007            NA -2.1848318
## tacs          -0.8158398      1.470842     -2.184832         NA
## [1] "###Condition by time z-scores###"
##                      tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                   NA    0.09861327    -0.1295825 0.3083629
## neurocaps_mcr  0.09861327            NA    -0.1923227 0.1838783
## neurocaps_ocr -0.12958247   -0.19232270            NA 0.3477431
## tacs           0.30836289    0.18387829     0.3477431        NA
## [1] "###Condition p-values###"
##                    tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                 NA  0.0170621969  0.1107349440 0.41459174
## neurocaps_mcr 0.0170622            NA  0.0004722765 0.14133389
## neurocaps_ocr 0.1107349  0.0004722765            NA 0.02890118
## tacs          0.4145917  0.1413338910  0.0289011789         NA
## [1] "###Condition by time p-values###"
##                    tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                 NA     0.9214453     0.8968968 0.7578062
## neurocaps_mcr 0.9214453            NA     0.8474894 0.8541089
## neurocaps_ocr 0.8968968     0.8474894            NA 0.7280331
## tacs          0.7578062     0.8541089     0.7280331        NA
## Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

##      ICC1  ICC1_lower ICC1_upper  ICC3   ICC3_lower ICC3_upper ICC3k ICC3k_lower ICC3k_upper               var
## 1   0.010 -0.39852452  0.4190756 0.010 -0.404808316  0.4216154 0.020 -1.36026200   0.5931497         -0.5 drug
## 2  -0.037 -0.43731361  0.3795720 0.000 -0.413247031  0.4132470 0.000 -1.40858948   0.5848192         -1.5 drug
## 3   0.000 -0.40701472  0.4106858 0.000 -0.413247031  0.4132470 0.000 -1.40858948   0.5848192          0.5 drug
## 4   0.439  0.03843318  0.7196487 0.439  0.030949557  0.7211326 0.610  0.06004088   0.8379745          1.5 drug
## 5   0.036 -0.37652310  0.4401861 0.036 -0.382935331  0.4426698 0.070 -1.24115137   0.6136814              drug
## 6  -0.012 -0.41685610  0.4007843 0.036 -0.383271168  0.4423533 0.069 -1.24291633   0.6133772      -0.5 neutral
## 7  -0.140 -0.51762725  0.2870705 0.064 -0.358693986  0.4649868 0.120 -1.11863596   0.6348000      -1.5 neutral
## 8   0.471  0.07975234  0.7390515 0.471  0.072303002  0.7404481 0.641  0.13485554   0.8508707       0.5 neutral
## 9   0.325 -0.09494472  0.6488050 0.359 -0.063215014  0.6727456 0.529 -0.13496163   0.8043609       1.5 neutral
## 10  0.355 -0.06111637  0.6680822 0.564  0.196853469  0.7926416 0.721  0.32895166   0.8843280           neutral
## 11  0.021 -0.38937000  0.4279682 0.033 -0.385178854  0.4405514 0.065 -1.25297855   0.6116427  -0.5DrugvNeutral
## 12 -0.012 -0.41669523  0.4009477 0.000 -0.413247031  0.4132470 0.000 -1.40858948   0.5848192  -1.5DrugvNeutral
## 13  0.188 -0.23692188  0.5559509 0.188 -0.243981312  0.5580785 0.317 -0.64543725   0.7163677   0.5DrugvNeutral
## 14  0.372 -0.04139961  0.6788837 0.420  0.007584464  0.7097232 0.591  0.01505475   0.8302200   1.5DrugvNeutral
## 15 -0.055 -0.45212917  0.3636542 0.000 -0.413247031  0.4132470 0.000 -1.40858948   0.5848192 DrugvNeutralEarly
## 16  0.116 -0.30521874  0.5029143 0.125 -0.303592532  0.5120978 0.223 -0.87188189   0.6773342  DrugvNeutralLate
## 17 -0.054 -0.45070979  0.3652001 0.000 -0.413247031  0.4132470 0.000 -1.40858948   0.5848192      DrugvNeutral
## Saving 7 x 5 in image

## [1] "###tDCS MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time +  (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 1320.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1302 -0.6149 -0.2205  0.4779  3.4277 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.1413   0.3759  
##  id           (Intercept) 0.1186   0.3444  
##  Residual                 0.6601   0.8125  
## Number of obs: 494, groups:  condition:id, 129; id, 65
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          2.09015    0.08225 119.66236  25.411  < 2e-16 ***
## conditiondrug       -0.24541    0.09890  64.12233  -2.481   0.0157 *  
## time                -0.22534    0.04734 373.25097  -4.760 2.78e-06 ***
## conditiondrug:time   0.09599    0.06602 370.24830   1.454   0.1468    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.606              
## time        -0.019  0.015       
## cndtndrg:tm  0.014 -0.018 -0.717
## [1] "###Neurocaps MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time +  (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 570.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8201 -0.5548 -0.2002  0.3057  4.1821 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.06786  0.2605  
##  id           (Intercept) 0.14683  0.3832  
##  Residual                 0.55775  0.7468  
## Number of obs: 228, groups:  condition:id, 58; id, 29
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          1.52463    0.11116  48.04003  13.716  < 2e-16 ***
## conditiondrug        0.12622    0.12038  27.94125   1.049    0.303    
## time                -0.29659    0.06221 168.51298  -4.768 4.01e-06 ***
## conditiondrug:time   0.08512    0.08824 168.78397   0.965    0.336    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.545              
## time        -0.003  0.003       
## cndtndrg:tm  0.002 -0.007 -0.705
## [1] "###Neurocaps OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time +  (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 481.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4627 -0.6814 -0.1119  0.5459  2.8380 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.02094  0.1447  
##  id           (Intercept) 0.19729  0.4442  
##  Residual                 0.76429  0.8742  
## Number of obs: 174, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          2.23722    0.13681  34.31196  16.352  < 2e-16 ***
## conditiondrug       -0.51833    0.13965  20.74077  -3.712  0.00131 ** 
## time                -0.44598    0.08432 128.93920  -5.289 5.12e-07 ***
## conditiondrug:time   0.11365    0.11926 128.92380   0.953  0.34237    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.510              
## time        -0.012  0.012       
## cndtndrg:tm  0.008 -0.016 -0.707
## [1] "###tACS OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time +  (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 391.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6610 -0.6255 -0.2142  0.4807  2.8207 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.04606  0.2146  
##  id           (Intercept) 0.23621  0.4860  
##  Residual                 0.40815  0.6389  
## Number of obs: 173, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          1.75887    0.13240  30.35512  13.284 3.47e-14 ***
## conditiondrug       -0.12052    0.11684  21.27910  -1.032  0.31389    
## time                -0.20351    0.06099 127.29143  -3.337  0.00111 ** 
## conditiondrug:time   0.06236    0.08679 127.67725   0.719  0.47372    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.439              
## time        -0.003  0.004       
## cndtndrg:tm  0.002 -0.008 -0.703

Box RT

plot_combination('box_rt', 'CravingRT', limits_y = c(0.5, 3.5), imaging = FALSE)
## [1] "###Condition z-scores###"
##                    tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                 NA     -0.789830     2.0756962  1.6805754
## neurocaps_mcr -0.789830            NA     2.4623103  2.1387617
## neurocaps_ocr  2.075696      2.462310            NA -0.6546852
## tacs           1.680575      2.138762    -0.6546852         NA
## [1] "###Condition by time z-scores###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA      2.534853      1.091190 -0.1697659
## neurocaps_mcr  2.5348532            NA     -1.416087 -2.4481921
## neurocaps_ocr  1.0911901     -1.416087            NA -1.1194922
## tacs          -0.1697659     -2.448192     -1.119492         NA
## [1] "###Condition p-values###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA    0.42962703    0.03792206 0.09284541
## neurocaps_mcr 0.42962703            NA    0.01380452 0.03245497
## neurocaps_ocr 0.03792206    0.01380452            NA 0.51267043
## tacs          0.09284541    0.03245497    0.51267043         NA
## [1] "###Condition by time p-values###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA    0.01124945     0.2751893 0.86519428
## neurocaps_mcr 0.01124945            NA     0.1567501 0.01435751
## neurocaps_ocr 0.27518926    0.15675014            NA 0.26293022
## tacs          0.86519428    0.01435751     0.2629302         NA
## Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

##      ICC1  ICC1_lower ICC1_upper  ICC3  ICC3_lower ICC3_upper ICC3k ICC3k_lower ICC3k_upper               var
## 1   0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000  -1.4085895   0.5848192         -0.5 drug
## 2   0.671  0.36328604  0.8480396 0.671  0.35676460  0.8489036 0.803   0.5259049   0.9182778         -1.5 drug
## 3   0.195 -0.23004587  0.5609553 0.292 -0.13806489  0.6291989 0.452  -0.3203603   0.7724028          0.5 drug
## 4   0.045 -0.36852838  0.4476398 0.049 -0.37193485  0.4529254 0.093  -1.1843830   0.6234668          1.5 drug
## 5   0.498  0.11444423  0.7545390 0.498  0.10704348  0.7558644 0.665   0.1933862   0.8609599              drug
## 6   0.276 -0.14787226  0.6165991 0.276 -0.15446338  0.6189687 0.433  -0.3653618   0.7646457      -0.5 neutral
## 7   0.075 -0.34282183  0.4708619 0.075 -0.34941686  0.4732592 0.139  -1.0741651   0.6424657      -1.5 neutral
## 8   0.810  0.60106781  0.9159504 0.810  0.59626054  0.9164458 0.895   0.7470717   0.9564015       0.5 neutral
## 9   0.249 -0.17599277  0.5983798 0.306 -0.12247901  0.6386818 0.469  -0.2791478   0.7795068       1.5 neutral
## 10  0.365 -0.04945140  0.6745100 0.365 -0.05692284  0.6761878 0.535  -0.1207173   0.8068163           neutral
## 11  0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000  -1.4085895   0.5848192  -0.5DrugvNeutral
## 12  0.715  0.43418603  0.8701040 0.715  0.42808612  0.8708513 0.834   0.5995242   0.9309680  -1.5DrugvNeutral
## 13 -0.114 -0.49824159  0.3108622 0.000 -0.41324703  0.4132470 0.000  -1.4085895   0.5848192   0.5DrugvNeutral
## 14  0.325 -0.09457070  0.6490235 0.457  0.05410178  0.7320805 0.627   0.1026500   0.8453192   1.5DrugvNeutral
## 15  0.224 -0.20099234  0.5814810 0.227 -0.20541124  0.5854176 0.370  -0.5170253   0.7385027 DrugvNeutralEarly
## 16  0.467  0.07401810  0.7364226 0.467  0.06656248  0.7378312 0.637   0.1248168   0.8491403  DrugvNeutralLate
## 17  0.349 -0.06786693  0.6643120 0.349 -0.07532117  0.6660319 0.517  -0.1629131   0.7995428      DrugvNeutral
## Saving 7 x 5 in image

## [1] "###tDCS MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time +  (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 1086.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8130 -0.5455 -0.1869  0.2686  3.4076 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.1706   0.4130  
##  id           (Intercept) 0.1763   0.4199  
##  Residual                 0.6704   0.8188  
## Number of obs: 396, groups:  condition:id, 121; id, 65
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          1.38757    0.09252  99.37271  14.998  < 2e-16 ***
## conditiondrug        0.67315    0.11522  56.46129   5.842 2.68e-07 ***
## time                 0.10055    0.04915 287.77919   2.045   0.0417 *  
## conditiondrug:time  -0.11214    0.07481 291.39762  -1.499   0.1350    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.539              
## time         0.078 -0.053       
## cndtndrg:tm -0.051  0.058 -0.656
## [1] "###Neurocaps MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time +  (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 428.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8776 -0.6184 -0.1551  0.2547  3.1725 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.08263  0.2875  
##  id           (Intercept) 0.14869  0.3856  
##  Residual                 0.71649  0.8465  
## Number of obs: 156, groups:  condition:id, 48; id, 26
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          1.15054    0.12925  30.42153   8.902 5.65e-10 ***
## conditiondrug        0.83342    0.16704  20.31452   4.989 6.74e-05 ***
## time                -0.05272    0.07641 102.21618  -0.690 0.491729    
## conditiondrug:time  -0.48453    0.12643 112.93184  -3.832 0.000209 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.498              
## time         0.029 -0.018       
## cndtndrg:tm -0.021 -0.043 -0.607
## [1] "###Neurocaps OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time +  (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 346
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9132 -0.4724 -0.1693  0.3446  3.3533 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.1120   0.3347  
##  id           (Intercept) 0.3614   0.6012  
##  Residual                 0.4676   0.6838  
## Number of obs: 139, groups:  condition:id, 40; id, 20
## 
## Fixed effects:
##                    Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)         1.52428    0.17328 25.16252   8.797 3.77e-09 ***
## conditiondrug       0.26377    0.16006 15.67478   1.648   0.1193    
## time                0.03842    0.07248 93.95191   0.530   0.5973    
## conditiondrug:time -0.25231    0.10443 97.08704  -2.416   0.0176 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.428              
## time         0.001  0.000       
## cndtndrg:tm  0.002  0.026 -0.693
## [1] "###tACS OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time +  (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 408.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6896 -0.5276 -0.1386  0.2052  4.5287 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.02208  0.1486  
##  id           (Intercept) 0.42820  0.6544  
##  Residual                 0.49879  0.7063  
## Number of obs: 165, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          1.27609    0.16272  25.57989   7.842 2.87e-08 ***
## conditiondrug        0.39445    0.11926  15.04231   3.307  0.00477 ** 
## time                -0.04465    0.06858 111.57509  -0.651  0.51634    
## conditiondrug:time  -0.09104    0.09924 114.56134  -0.917  0.36086    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.360              
## time         0.008 -0.007       
## cndtndrg:tm -0.005  0.016 -0.689
extract_slopes_beh <- function(measure, this_label, idps, prefix = 'dcr_'){
  #same as above, but for behavioral data
  #prefix will be dcr_tpre_ for tdcs/tacs and dcr_ for neurocaps 
  #measure will be box_rt, craving_rt, or craving
  #measure <- 'craving'
  #prefix <- 'dcr_tpre_'
  #this_label <- 'Craving'
  #idps <- idps_tdcs
  
  
  #idps_tacs$dcr_tpre_response_craving_[0-7]
  #idp_tdcs
  #dcr_tpre_response_box_rt_0
  #dcr_tpre_response_craving_rt_0
  #idps_neurocaps_ocr$dcr_response_box_craving[0-7]

  column_suffixes <- 0:7
  
  these_cols <- c('id', paste0(prefix, 'response_', measure, '_', column_suffixes))
  library(reshape2)
  one_dataset <- idps[, these_cols]
  long_data <- melt(one_dataset, id.vars = c('id'))
  long_data$variable <- as.character(long_data$variable)
  long_data$condition <- NA
  long_data$number <- substr(long_data$variable, nchar(long_data$variable), nchar(long_data$variable))
  long_data$condition[long_data$number %in% c('0', '2', '4', '6')] <- 'neutral'
  long_data$condition[long_data$number %in% c('1', '3', '5', '7')] <- 'drug'
  #put neutral in the intercept
  long_data$condition <- factor(long_data$condition, levels = c('neutral', 'drug'))
  #time = block number, just like for the imaging variables
  long_data$time <- NA
  long_data$time[long_data$number %in% c('0', '1')] <- 1
  long_data$time[long_data$number %in% c('2', '3')] <- 2
  long_data$time[long_data$number %in% c('4', '5')] <- 3
  long_data$time[long_data$number %in% c('6', '7')] <- 4

  #mean center on time
  long_data$time <- long_data$time - mean(long_data$time)

  library(lme4)
  library(lmerTest)
  library(ggplot2)
  library(sjstats) #for icc of mixed effects models
  
  
  
  all_results <- NULL
  for (this_id in unique(long_data$id)){  
    #this_id <- 'AR316'
    this_sub_data <- long_data[long_data$id == this_id,]
    neutral_string <- paste0(this_label, '_slopeNeutral')
    meth_string <- paste0(this_label, '_slopeMeth')
    delta_string <- paste0(this_label, '_slopeMethMinusNeutral')
    if (sum(!(is.na(this_sub_data$value[this_sub_data$condition == 'neutral']))) > 1 &
        sum(!(is.na(this_sub_data$value[this_sub_data$condition == 'drug']))) > 1){
      #subject made at least two ratings per condition
      this_lm <- lm(paste('value ~ condition * time'), data = this_sub_data)
      
      these_results <- data.frame(list(id = this_id, neutral_string = this_lm$coefficients['time'],
                                meth_string = this_lm$coefficients['time'] + this_lm$coefficients['conditiondrug:time'],
                                delta_string = this_lm$coefficients['conditiondrug:time']))
      names(these_results) <- c('id', neutral_string, meth_string, delta_string)
      
    } else{
      #not enough ratings to estimate two slopes
      these_results <- data.frame(list(id = this_id, neutral_string = NA,
                                meth_string = NA,
                                delta_string = NA))
      names(these_results) <- c('id', neutral_string, meth_string, delta_string)

        }
    all_results <- rbind(all_results, these_results)
  }
  

  return(all_results)
   
}

Corrplot including individual slopes and self-report

extract_one_slopes <- function( roi = '1', this_label = 'VMPFC', prefix = 'dcr_tpre_', idps = idps_tdcs){
  column_prefixes <- c('stats_tdcsprelim_drug.r11.0.coef_mean_',
           'stats_tdcsprelim_drug.r12.0.coef_mean_',
           'stats_tdcsprelim_drug.r13.0.coef_mean_',
           'stats_tdcsprelim_drug.r14.0.coef_mean_',
           'stats_tdcsprelim_neutral.r11.0.coef_mean_',
           'stats_tdcsprelim_neutral.r12.0.coef_mean_',
           'stats_tdcsprelim_neutral.r13.0.coef_mean_',
           'stats_tdcsprelim_neutral.r14.0.coef_mean_')
  this_roi <- c('id', 'motion', paste0(prefix, column_prefixes, roi))
  one_dataset <- idps[, this_roi]
  long_data <- melt(one_dataset, id.vars = c('id', 'motion'))
  long_data$condition <- NA
  long_data$condition[grepl('neutral', long_data$variable)] <- 'neutral'
  long_data$condition[grepl('drug', long_data$variable)] <- 'drug'
  #put neutral in the intercept
  long_data$condition <- factor(long_data$condition, levels = c('neutral', 'drug'))
  long_data$time <- NA
  long_data$time[grepl('r11', long_data$variable)] <- 1
  long_data$time[grepl('r12', long_data$variable)] <- 2
  long_data$time[grepl('r13', long_data$variable)] <- 3
  long_data$time[grepl('r14', long_data$variable)] <- 4
  #mean center on time
  long_data$time <- long_data$time - mean(long_data$time)
  
  #this_lme <- lmer(paste('value ~ condition * time + (1|id)'), data = long_data)
  #for checking for NA's, looks like it's all good now
  #print(long_data)
  all_results <- NULL
  for (this_id in unique(long_data$id)){  
    #this_id <- 'AR316'
    this_lm <- lm(paste('value ~ condition * time + motion'), data = long_data[long_data$id == this_id,])
    neutral_string <- paste0(this_label, '_slopeNeutral')
    meth_string <- paste0(this_label, '_slopeMeth')
    delta_string <- paste0(this_label, '_slopeMethMinusNeutral')
    these_results <- data.frame(list(id = this_id, neutral_string = this_lm$coefficients['time'],
                                meth_string = this_lm$coefficients['time'] + this_lm$coefficients['conditiondrug:time'],
                                delta_string = this_lm$coefficients['conditiondrug:time']))
    names(these_results) <- c('id', neutral_string, meth_string, delta_string)
    all_results <- rbind(all_results, these_results)
  }
  
  

  return(all_results)
}


all_slopes <- extract_one_slopes('1', 'VMPFC', prefix = 'dcr_tpre_', idps = idps_tdcs)
all_slopes <- merge(all_slopes, extract_one_slopes('3', 'LSTG', prefix = 'dcr_tpre_', idps = idps_tdcs))
all_slopes <- merge(all_slopes, extract_one_slopes('4', 'RSTG', prefix = 'dcr_tpre_', idps = idps_tdcs))
all_slopes <- merge(all_slopes, extract_one_slopes('5', 'LVSTRI', prefix = 'dcr_tpre_', idps = idps_tdcs))
all_slopes <- merge(all_slopes, extract_one_slopes('6', 'RVSTRI', prefix = 'dcr_tpre_', idps = idps_tdcs))
all_slopes <- merge(all_slopes, extract_one_slopes('8', 'RAmy', prefix = 'dcr_tpre_', idps = idps_tdcs))






vas_data <- read.csv('../paper-dcr-temporaldynamics/data/MethVASData.csv')
vas_data <- vas_data[, c('record_id', 'redcap_event_name', 'mcs_vas', 'mcs_vas_2')]

names(vas_data) <- c('id', 'visit', 'craving', 'control')

library(reshape2)
vas_data_1 <- vas_data[vas_data$visit == 'before_pre_fmri_arm_1', c('id', 'craving', 'control')]
vas_data_2 <- vas_data[vas_data$visit == 'after_pre_fmri_arm_1', c('id', 'craving', 'control')]

vas_data_wide <- merge(vas_data_1, vas_data_2, by = 'id')

names(vas_data_wide) <- c('id', 'craving_pre', 'control_pre', 'craving_post', 'control_post')
vas_data_wide$craving_delta <- vas_data_wide$craving_post - vas_data_wide$craving_pre
vas_data_wide$control_delta <- vas_data_wide$control_post - vas_data_wide$control_pre

merged_data <- merge(all_slopes, vas_data_wide, all.x = TRUE)



last_use_data <- read.csv('../paper-dcr-temporaldynamics/data/ADUQ_2019-06-12_1407.csv')
merged_data <- merge(merged_data, last_use_data[, c('id', 'aduq_20a')], all.x = TRUE)
names(merged_data)[names(merged_data) == 'aduq_20a'] <- 'DaysSinceLastUse'




tableone_data <- read.csv('../paper-dcr-temporaldynamics/data/Table1Database66_ver2.csv')
merged_data <- merge(merged_data, tableone_data, all.x = TRUE)

idps_tdcs$methvsneutral_craving_selfreport_insidescanner <- (idps_tdcs$dcr_tpre_response_craving_1 + idps_tdcs$dcr_tpre_response_craving_3 +
  idps_tdcs$dcr_tpre_response_craving_5 + idps_tdcs$dcr_tpre_response_craving_7 - 
    idps_tdcs$dcr_tpre_response_craving_0 - idps_tdcs$dcr_tpre_response_craving_2 - idps_tdcs$dcr_tpre_response_craving_4 - idps_tdcs$dcr_tpre_response_craving_6) / 4


merged_data <- merge(merged_data, extract_slopes_beh('craving', 'Craving', idps_tdcs, prefix = 'dcr_tpre_'), all.x = TRUE)

merged_data <- merge(merged_data, idps_tdcs[, c('id', 'methvsneutral_craving_selfreport_insidescanner')])
names(merged_data)
##  [1] "id"                                                "VMPFC_slopeNeutral"                                "VMPFC_slopeMeth"                                   "VMPFC_slopeMethMinusNeutral"                       "LSTG_slopeNeutral"                                
##  [6] "LSTG_slopeMeth"                                    "LSTG_slopeMethMinusNeutral"                        "RSTG_slopeNeutral"                                 "RSTG_slopeMeth"                                    "RSTG_slopeMethMinusNeutral"                       
## [11] "LVSTRI_slopeNeutral"                               "LVSTRI_slopeMeth"                                  "LVSTRI_slopeMethMinusNeutral"                      "RVSTRI_slopeNeutral"                               "RVSTRI_slopeMeth"                                 
## [16] "RVSTRI_slopeMethMinusNeutral"                      "RAmy_slopeNeutral"                                 "RAmy_slopeMeth"                                    "RAmy_slopeMethMinusNeutral"                        "craving_pre"                                      
## [21] "control_pre"                                       "craving_post"                                      "control_post"                                      "craving_delta"                                     "control_delta"                                    
## [26] "DaysSinceLastUse"                                  "Age"                                               "Education"                                         "BMI"                                               "Age.of.Meth.use.onset..years.old."                
## [31] "Duration.of.Meth.use.at.least.once.a.week..years." "Cost.of.Meth..dollar.per.month."                   "Dose.of.Meth..gram.per.day."                       "History.of.Meth.injection.b..n....."               "Days.of.Meth"                                     
## [36] "Days.of.Alcohol"                                   "Days.of.Alcohol.intoxication"                      "Days.of.Heroin"                                    "Days.of.Methadone"                                 "Days.of.Other.opioids"                            
## [41] "Days.of.Barbiturate"                               "Days.of.Sedative"                                  "Days.of.Cocaine"                                   "Days.of.Cannabis"                                  "Days.of.Hallucinogens"                            
## [46] "Days.of.Inhalants"                                 "Duration.of.current.abstinence..days."             "Meth.Cue.Reactivity.Screening.score..0.100."       "Meth.Withdrawal.Scale.score..0.24."                "Motion..Rest.fMRI..Pre.tDCS"                      
## [51] "Motion..Task.fMRI..Pre.tDCS"                       "Motion..Rest.fMRI..Post.tDCS"                      "Motion..Task.fMRI..Post.tDCS"                      "Craving_slopeNeutral"                              "Craving_slopeMeth"                                
## [56] "Craving_slopeMethMinusNeutral"                     "methvsneutral_craving_selfreport_insidescanner"
names(merged_data)[names(merged_data) == 'Meth.Cue.Reactivity.Screening.score..0.100.'] <- 'BaselineCueReactivity'
names(merged_data)[names(merged_data) == 'Duration.of.Meth.use.at.least.once.a.week..years.'] <- 'MethUseDuration'
names(merged_data)[names(merged_data) == 'Cost.of.Meth..dollar.per.month.'] <- 'MethCost'

names(merged_data)[names(merged_data) == 'craving_post'] <- 'Craving_post'
names(merged_data)[names(merged_data) == 'craving_pre'] <- 'Craving_pre'
names(merged_data)[names(merged_data) == 'craving_delta'] <- 'Craving_delta'

to_plot <- c('Age', 'MethUseDuration', 'MethCost',"DaysSinceLastUse",
            'BaselineCueReactivity',  "Craving_pre", "Craving_post", "Craving_delta",
            "VMPFC_slopeNeutral", "VMPFC_slopeMeth",  "VMPFC_slopeMethMinusNeutral",
            "LSTG_slopeNeutral", "LSTG_slopeMeth","LSTG_slopeMethMinusNeutral",
            "RSTG_slopeNeutral", "RSTG_slopeMeth", "RSTG_slopeMethMinusNeutral",
            "LVSTRI_slopeNeutral", "LVSTRI_slopeMeth", "LVSTRI_slopeMethMinusNeutral",
            "RVSTRI_slopeNeutral", "RVSTRI_slopeMeth", "RVSTRI_slopeMethMinusNeutral",
            "RAmy_slopeNeutral", "RAmy_slopeMeth", "RAmy_slopeMethMinusNeutral")

to_plot <- c('Age', 'MethUseDuration', 'MethCost',"DaysSinceLastUse",
            "Craving_pre", "Craving_post", "Craving_delta", 
            'methvsneutral_craving_selfreport_insidescanner',
            'Craving_slopeMethMinusNeutral',
            "VMPFC_slopeMethMinusNeutral",
            "LSTG_slopeMethMinusNeutral",
            "RSTG_slopeMethMinusNeutral",
            "LVSTRI_slopeMethMinusNeutral",
            "RVSTRI_slopeMethMinusNeutral",
            "RAmy_slopeMethMinusNeutral")

library(corrplot)
## corrplot 0.84 loaded
this_matrix <- cor(merged_data[, to_plot], use = 'pairwise.complete.obs')
pvals_raw <- cor.mtest(merged_data[, to_plot])

#get FDR corrected p-values to use for plotting
pvals_fdr <- pvals_raw$p
pvals_fdr[upper.tri(pvals_fdr)] <- p.adjust(pvals_raw$p[upper.tri(pvals_raw$p)])
pvals_fdr[lower.tri(pvals_fdr)] <- p.adjust(pvals_raw$p[lower.tri(pvals_raw$p)])


col2 <- colorRampPalette(rev(c("#67001F", "#B2182B", "#D6604D", "#F4A582",
                           "#FDDBC7", "#FFFFFF", "#D1E5F0", "#92C5DE",
                           "#4393C3", "#2166AC", "#053061")))


png('corrplot_slopes_sxs.png', width = 1000, height = 1000)
corrplot.mixed(this_matrix, upper = 'ellipse', lower = 'number', tl.pos = 'lt', tl.cex = 1,#, lower.col = "black",
  upper.col = col2(50), lower.col = col2(50), p.mat = pvals_fdr, sig.level = 0.05)
dev.off()
## png 
##   2